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ABSTRACT 

A  formula  for  the  tidal  dissipation  rate  in  a  spherical  body  is  derived  from  first  principles  to  correct  some 
mathematical  inaccuracies  found  in  the  literature.  The  development  is  combined  with  the  Darwin-Kaula  formalism 
for  tides.  Our  intermediate  results  are  compared  with  those  by  Zschau  and  Platzman.  When  restricted  to  the  special 
case  of  an  incompressible  spherical  planet  spinning  synchronously  without  libration,  our  final  formula  can  be 
compared  with  the  commonly  used  expression  from  Peale  &  Cassen.  However,  the  two  turn  out  to  differ,  as  in  our 
expression  the  contributions  from  all  Fourier  modes  are  positive-definite,  which  is  not  the  case  with  the  formula 
from  Peale  &  Cassen.  Examples  of  the  application  of  our  expression  for  the  tidal  damping  rate  are  provided  in  the 
work  by  Makarov  &  Efroimsky  (Paper  II)  published  back  to  back  with  the  current  paper. 
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1.  MOTIVATION  AND  PLAN 

The  tidal  heating  of  planets  and  moons  has  long  been  a 
key  area  of  planetary  science.  Accurate  investigation  into  this 
process  requires  numerical  integration  of  dissipation  over  layers 
of  the  perturbed  body.  At  the  same  time,  it  is  common  to 
infer  qualitative  conclusions  from  approximations  based  on 
modeling  the  body  with  a  homogeneous  sphere  of  a  certain 
rheology.  However,  the  simplistic  nature  of  the  approach  limits 
the  precision  of  the  ensuing  conclusions.  For  example,  the 
presence  of  a  sizable  molten  core,  such  as  that  in  Mercury,  may 
increase  the  damping  rate  compared  to  a  homogeneous  body. 
Still,  estimates  obtained  with  our  simplified  homogeneous- 
sphere  model  should  be  accurate  within  a  factor  of  several — thus 
(1)  serving  as  a  useful  guidance  for  solar  system  bodies,  and  (2) 
being  completely  legitimate  for  exoplanets,  as  our  knowledge 
of  their  structure  is  speculative  at  best. 

In  our  paper,  we  derive  from  the  first  principles  a  formula  for 
the  tidal  heating  rate  in  a  tidally  perturbed  homogeneous  sphere. 
We  compare  our  result  with  the  formulae  used  in  the  literature 
and  point  out  the  differences. 

In  Sections  3-5,  we  present  an  accurate  re-examination  of 
the  standard  integral  expression  for  the  damping  rate  in  a 
homogeneous  incompressible  sphere  subject  to  tides.  The  check 
is  necessary  because  in  previous  studies  the  expression  was 
derived  in  an  ad  hoc  manner,  sometimes  with  demonstrable 
mathematical  inaccuracies.  The  conventional  derivation  begins 
with  the  general  formula  for  the  power 

P  —  J  pEv  ■  VV'E  re¬ 
written  in  the  Eulerian  description  (i.e.,  via  coordinates  associ¬ 
ated  with  a  deformed  body).  Its  time  average  is  then  cast  into 
the  form  of 

1  00  r 

<p>  =  4 ^D2i+1>/WTO<iS 


which  is  in  the  Lagrangian  language  (an  integral  over  an 
undeformed  body).  In  the  former  equation,  pe  is  the  Eulerian 
density,  v  is  the  Eulerian  velocity,  V'E  is  the  Eulerian  perturbation 
of  the  potential  (perturbation  assembled  of  the  tide-raising 
potential  and  the  resulting  additional  tidal  potential  of  the 
deformed  body),  and  r  is  a  perturbed  position  in  the  body  frame. 
In  the  latter  equation,  Wi  and  U )  are  the  degree-/  components  of 
the  tide-raising  and  additional  tidal  potentials,  G  is  the  Newton 
gravity  constant,  R  is  the  radius  of  the  planet,  and  dS  is  an 
element  of  the  undeformed  surface  of  the  sphere. 

The  transition  from  the  former  formula  to  the  latter  requires 
the  use  of  the  boundary  conditions  on  the  free  surface.  At  that 
point,  integration  is  already  carried  out  within  the  Lagrangian 
description  (over  an  undeformed  surface),  but  the  boundary 
conditions  are  nonetheless  imposed  on  the  Eulerian  potential 
and  its  gradient.  (The  boundary  conditions  are  much  simpler  in 
the  Eulerian  form.)  This  mixed  treatment  requires  attention,  and 
its  employment  by  the  early  authors  (Zschau  1978;  Platzman 
1984)  contained  inaccuracies.  However,  none  of  those  turned 
out  to  be  critical,  and  the  above  expression  for  the  average 
power  (P)  is  correct  for  small  deformations. 

In  Section  6,  we  explore  the  standard  way  of  casting  the  above 
integral  into  a  spectral  sum  over  the  tidal  Fourier  modes  a>.  It  is 
commonly  assumed  that  the  result  should  read  as  in  Platzman 
(1984): 

(P)  =  4Jgr^2(21  +  1)  ^  W,2(w) £,(&>) sin e,(w). 


Here  kj(co)  and  e/(co)  are  the  Love  number  and  phase  lag 
corresponding  to  the  Fourier  mode  o>  =  coimpq ,  with  Impq  being 
the  four  integers  wherewith  the  Fourier  modes  are  numbered  in 
the  Darwin-Kaula  theory  of  tides  (see  Efroimsky  &  Makarov 
2013  and  references  therein).  However,  an  accurate  investigation 
demonstrates  that  the  spectral  sum  differs  from  the  above.  The 
difference  originates  for  two  reasons.  One  is  the  degeneracy, 
that  is,  the  fact  that  several  different  Fourier  modes  ooimpq  share  a 
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numerical  value  <y,  so  the  structure  of  the  sum  is  more  complex. 1 
The  second  reason  is  that  the  modes  can  be  of  either  sign, 
not  necessarily  positive.  So  the  resulting  power  will  contain 
seemingly  strange  terms  with  Wi(a>)  Wi(—co)ki(a>)  sine/(&)). 

These  difficulties  were  noticed  and  analyzed  by  Peale  & 
Cassen  (1978),  but  their  result  needs  correction  too.  Some 
terms  in  their  spectral  sum  (Equation  (31)  in  Peale  &  Cassen) 
are  not  positive  definite,  whence  an  underestimate  of  the  heat 
production  may  result.2 

The  calculation  of  the  power  production,  developed  by  Peale 
&  Cassen  (1978),  implies  averaging  not  only  over  the  tidal 
period,  but  also  over  the  apsidal  period.  This  can  be  observed 
from  the  formulae  (20)  and  (21)  in  their  work.  In  our  paper, 
however,  we  consider  two  separate  cases:  those  with  and  without 
apsidal  precession.  In  the  first  case,  the  period  of  the  apsidal 
precession  is  shorter  than  the  typical  time  of  relaxation  in 
the  mantle  (which  may  be  identified  with  the  Maxwell  time). 
The  argument  of  the  pericenter  of  the  perturber,  of,  cannot  be 
treated  as  constant,  so  the  formula  for  the  mean  power  should  be 
averaged  not  only  over  the  tidal  period,  but  also  over  the  period 
of  the  pericenter  motion.  (We  assume  this  motion  steady.)  In  the 
second  case,  the  evolution  of  the  line  of  apsides  is  slow,  with 
its  period  being  longer  than  the  Maxwell  time.  The  argument  of 
the  pericenter  should  be  regarded  as  a  constant.  Accordingly,  in 
the  latter  case  the  tidal  dissipation  formula  is  more  complicated 
because  it  includes  explicit  dependence  of  Fourier  terms  on  the 
argument  of  pericenter. 

This  paper  is  a  condensed  version  of  a  comprehensive  preprint 
(Efroimsky  &  Makarov  2014),  to  which  we  refer  for  technical 
details  and  derivations. 

In  a  subsequent  work,  Makarov  &  Efroimsky  (2014,  hereafter 
Paper  II),  published  back  to  back  with  this  one,  we  apply  our 
results  in  three  case  studies:  Io,  Mercury,  and  Kepler-  10b.  In  that 
paper  we,  among  other  things,  hypothesize  that  the  tidal  heating 
rate  at  spin-orbit  resonances  is  greatly  influenced  by  libration 
and,  therefore,  by  the  triaxiality  of  the  tidally  perturbed  body. 


Things  are  complicated  even  further  because  different  mecha¬ 
nisms  of  friction  become  leading  over  different  frequency  bands, 
thus  the  tidal  response  cannot  be  described  by  one  simple  dissi¬ 
pation  model  (Efroimsky  2012a,  2012b). 

2.1.  Generalities 

The  development  of  the  mathematical  theory  of  bodily  tides 
was  started  by  Darwin  (1879)  who  derived  a  partial  sum  of 
the  Fourier  expansion  of  the  additional  potential  of  a  tidally 
perturbed  sphere.  A  decisive  contribution  into  this  theory  was 
offered  almost  a  century  later  by  Kaula  ( 1 964),  who  wrote  down 
a  complete  series.  In  a  previous  paper  (Efroimsky  &  Makarov 
2013),  we  provided  a  detailed  presentation  of  the  Darwin-Kaula 
expansion,  and  explained  how  tidal  friction  and  lagging  are 
built  into  it.  We  compared  the  Darwin-Kaula  theory  with  the 
one  by  MacDonald  (1964),  and  demonstrated  that  the  former 
theory  is  superior  to  the  latter  because  it  can,  in  principle,  be 
combined  with  an  arbitrary  rheology.  Referring  the  reader  to  the 
aforementioned  literature  for  details,  we  present  several  central 
formulae  that  will  be  necessary. 

An  external  body  of  mass  M* ,  located  in  r  *  =  (r* ,  X* ,  (jf  ), 
generates  the  following  disturbing  potential  in  a  point  R  — 
(R,  <p,  X)  on  the  surface  of  a  sphere  of  radius  R  <  r*: 


W(R,  r*)  =  W,(R,r *) 

1=2 

GM 


GM* 


E(£ 


1=2 


P/icos  y) 


R 


E  S  E 


i  i 


(l  —  m)\ 
(/  +  m)\ 


(2  —  S0m)P/„ 


1=2  v  '  m= 0 

x  (sin0).P/m(sin0*)cos  777(.A.  —  X*).  (1) 


2.  THE  DARWIN-KAULA  FORMALISM  IN  BRIEF 

Describing  linear  bodily  tides  consists  of  two  steps.  First,  it  is 
necessary  to  Fourier-expand  both  the  tide-raising  potential  and 
the  induced  additional  potential  of  the  tidally  perturbed  body. 
Second,  it  is  necessary  to  link  each  Fourier  component  of  the 
additional  tidal  potential  to  an  appropriate  Fourier  component  of 
the  tide-raising  potential.  This  means  establishing  the  phase  lag 
and  the  ratio  of  magnitudes  called  the  dynamical  Love  number. 

Due  to  interplay  of  rheology  and  self-gravitation,  the  phase 
lags  and  Love  numbers  have  nontrivial  frequency  dependencies. 


1  When  calculating  Wp  one  has  first  to  group  together  and  sum  all  the  terms 
corresponding  to  a  particular  value  of  co.  Each  sum  should  be  squared  and 
averaged,  and  only  after  that  should  the  final  summation  over  the  distinct 
values  of  co  be  earned  out.  In  the  original  expression  for  the  average  power, 

22^(2/  +  1 )  {co/2)  W2(o>)  ki(co)  sin  c/(oj)  ,  the  W2{co)  term  should 
be  replaced  with  the  squared  sum  of  all  the  harmonics  of  W  that  correspond  to 
a  particular  value  of  co. 

2  In  the  expression  for  ( P) ,  an  input  from  each  value  of  coimpq  must  be 
non-negative.  This  can  be  observed  from  the  fact  that  the  mode  co  =  co\mpq  and 
the  corresponding  phase  lag  ci(co)  =  co  A//(7v)  are  always  of  the  same  sign  (the 
time  lag  A tj(co)  being  positive  definite  due  to  causality).  Thus  the  product 
co€i{co)  =  coimpq  ei(cojmpq)  in  the  spectral  sum  can  always  be  rewritten  as 
\<vimpq\/  Qlmpq,  with  the  tidal  quality  factor  being  defined  via 

1  / Qimpq  =  |  sin  cj((c)jmpq)\.  In  their  spectral  sum,  Peale  &  Cassen  (1978, 
Equation  (31))  have  just  a>impq/ Qimpq,  and  not  \(oimpq\/ Qimpq.  As  a  result, 
some  terms  come  out  negative  and  the  heat  production  intensity  may  be 
underestimated. 


Here  G  denotes  Newton’s  gravity  constant,  0  is  the  latitude 
reckoned  from  the  spherical  body’s  equator,  X  is  the  longitude 
measured  from  a  fixed  meridian,  and  y  is  the  angular  separation 
between  the  vectors  r  *  and  R  pointing  from  the  perturbed 
body’s  center.  The  definition  and  normalization  of  the  Legendre 
polynomials  P/( cos  y )  and  the  associated  Legendre  polynomials 
P/m(sin  0)  are  given  in  Appendix  A. 

While  in  the  above  formula  the  location  of  the  perturber  on 
its  trajectory  is  expressed  through  the  spherical  coordinates 
r*  —  (r* ,  X* ,  0*),  a  trigonometric  transformation  (developed 
by  Kaula  1961)  enables  one  to  switch  to  the  perturber’s  orbital 
elements  r*  =  (a*,  e* ,  i*.  LI* ,  of,  M*).  Thus,  the  disturbing 
potential  is  expressed  as 


W(R,  r*)=Y  W,mpq 

Impq 


GM * 


1=2 


E(£)E 


i  i 
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where  9*  is  the  rotation  angle  of  the  tidally  perturbed  body,3 
while  Fimp(i*)  and  Gipq(e*)  are  the  inclination  functions  and 
the  eccentricity  polynomials,  respectively.  The  auxiliary  linear 
combinations  v*  are  defined  by 

V*lmpq  =  {l  -  2 p)  «*+(/-  2p  +  q)M*  +  m  Q *.  (3) 

Conventionally,  the  letters  denoting  the  elements  of  the  per- 
turber  are  accompanied  with  asterisks:  a*,  e* ,  i* ,  Q*,  of ,  M*. 
Following  Kaula  (1964),  the  sidereal  angle  also  acquires  an  as¬ 
terisk  when  it  appears  in  a  combination  <npq  ~  m  With  the 
perturber’s  elements. 

The  angle  9,  however,  does  not  acquire  an  asterisk  when  it 
appears  in  a  linear  combination  Vimpq  —  m  9  with  the  orbital 
elements  of  a  test  body  subject  to  the  additional  tidal  potential 
of  the  perturbed  body.  This  strange  nomenclature  introduced 
by  Kaula  (1964) — two  different  notations  for  one  angle — turns 
out  to  be  helpful  and  convenient  in  the  calculation  of  the 
back-reaction  experienced  by  the  perturber.  For  comprehensive 
explanation  of  this  obscure  point,  see  Section  5  in  Efroimsky  & 
Makarov  (2013). 

Over  timescales  shorter  than  the  apsidal-motion  period, 
the  expression  in  round  brackets  in  the  formula  (2)  can  be 
linearized  as 

V*mpq  ~  =  Mlmpq  V  “  fo) 

-ml  +  vlpq(t0)  -  mO*(to),  (4) 

where  the  following  quantities  act  as  the  Fourier  tidal  modes: 

coimpq  =  V*[mpq  -  m  9*  =  (/  -  2 p)  cb*  +  (/  -  2p  +  q)  M  * 

+  m(Q*  -  9*),  (5) 

M.  *  being  the  perturber’s  “anomalistic”  mean  motion  (see 
Section  2.3  below),  and  to  being  the  time  of  pericenter  passage. 
(As  ever,  we  set  M.  *  =  0  in  the  pericenter.)  The  modes  u>impq 
can  assume  either  sign,  but  the  physical  forcing  frequencies  are 
positive  definite: 

Xlmpq  —  0)jinpq  | .  (6) 

2.2.  Simplifying  the  Notation:  Less  Asterisks 

In  the  preceding  subsection,  we  obeyed  the  convention  by 
Kaula  (1964)  and  marked  with  asterisk  the  orbital  elements 
of  the  tide-raising  body.  Kaula  introduced  this  notation  because 
within  his  model  he  also  considered  another  exterior  body,  which 
was  disturbed  by  the  tides  generated  on  the  planet  by  the  tide¬ 
raising  body.  This  exterior  body’s  elements  were  denoted  by  the 
same  letters,  but  without  an  asterisk. 

When  the  two  outer  bodies  coincide,  the  asterisks  may  be 
dropped,  except  on  two  occasions.  The  first  is  writing  the 
masses — while  the  mass  of  the  planet  is  denoted  with  M, 
the  mass  of  the  perturber  (the  star)  will  be  written  as  M*.  The 
other  occasion  requires  writing  the  additional  tidal  potential  of 
the  perturbed  body — the  potential  will  have  a  value  U(r,r*) 
in  a  point  r,  provided  the  perturber  (the  star)  is  located  in 
an  exterior  point  r  *  (both  vectors  being  planetocentric).  The 
planet’s  rotation  rate  9,  as  well  as  the  orbital  elements  of  the 
star  as  seen  from  the  planet,  will  hereafter  be  written  without 
asterisks. 

The  most  important  notations  employed  in  this  paper  are 
collected  in  Table  1 . 

3  When  the  equinoctial  precession  may  be  neglected,  8*  may  be  regarded  as 
the  sidereal  angle. 


Table  1 

Symbol  Key 


Notation 

Description 

r* 

The  position  of  the  star  relative  to  the  center  of  the 
planet 

r* 

The  star-planet  distance 

0* 

The  declination  of  the  star  relative  to  the  equator  of 
the  planet 

X* 

The  right  ascension  of  the  star  relative  to  a  fixed 
meridian  on  the  planet 

R 

A  point  on  the  surface  of  the  planet 

r 

A  point  outside  the  planet,  located  above  the  surface 
point  R 

R 

The  radius  of  the  planet 

r 

Distance  from  the  center  of  the  planet  to  an  exterior 
point  r 

0 

The  latitude  of  the  point  R  on  the  surface  of  the 
planet 

(also  the  declination  of  an  exterior  point  r  located 
above  the  surface  point  R) 

X 

The  longitude  of  a  point  R  on  the  surface  of  the 
planet 

(also  the  right  ascension  of  an  exterior  point  r 
located  above  the  surface  point  R) 

W(R.r*) 

Tide-raising  potential  at  a  surface  point  R  of  the 
planet 

W/(R.  r  *) 

The  /-degree  part  of  the  tide-raising  potential  at  a 
surface  point  R  of  the  planet 

U (r,  r *) 

Additional  tidal  potential  in  a  point  r  outside  the 
planet 

U,(r ,  r *) 

The  /-degree  part  of  the  additional  tidal  potential  in 
a  point  r  outside  the  planet 

V° 

The  constant-in-time  spherically  symmetrical 
potential  of  an  undeformed  planet 

V' 

The  total  perturbation  of  the  potential  of  the  planet 
(V'  =  W+  U ) 

V 

The  overall  potential  of  the  planet 

tv  =  v°+  v  =  y°+  w+  u) 

to) Impq 

The  Fourier  modes  of  the  tide 

Xlmpq 

The  physical  frequencies  of  the  tidal  stresses  and 
strains  (ximpq  =  \to)impq\) 

X 

A  shorter  notation  for  Xlmpq 

M 

The  mass  of  the  planet 

M* 

The  mass  of  the  star 

a 

The  semimajor  axis 

e 

The  eccentricity 

i 

The  obliquity  (the  inclination  of  the  star  as  seen 
from  the  planet) 

CO 

The  argument  of  the  pericenter  of  the  star  as  seen 
from  the  planet 

CO 

When  there  is  no  risk  of  confusion,  co  is  also  used  as 
a  short  notation  for  coimpq 

Q 

The  longitude  of  the  node  of  the  star  as  seen  from 
the  planet 

M 

The  mean  anomaly  of  the  star  as  seen  from  the 
planet 

n 

The  mean  motion 

G 

Newton’s  gravitational  constant 

8 

The  rotation  angle  of  the  planet 

8 

The  spin  rate  of  the  planet 

ki.hi 

The  degree-/  static  Love  numbers  of  the  planet 

(.to^lmpq  )  i  hl(to)lmpq ) 

The  degree-/  dynamical  Love  numbers  of  the  planet 

£/  (.to)  Impq) 

The  degree-/  tidal  phase  lag 

A  tl(to)impq) 

The  degree-/  tidal  time  lag 

Ql(to)[mpq) 

The  degree-/  tidal  quality  factor  defined  as 

Ql(to)lmpq )  —  1/ sin  \€l(to)impq)\ 

Qlmpq 

The  Q  factor,  in  the  notation  of  Kaula  (1964)  and 
Peale  &  Cassen  (1978) 

F Imp  (i ) 

The  inclination  functions 

3 


The  Astrophysical  Journal,  795:6  (19pp),  2014  November  1 


Efroimsky  &  Makarov 


Table  1 

(Continued) 


Notation 

Description 

Glpq  (^) 

The  eccentricity  polynomials 

u 

The  mean  rigidity  of  the  mantle 

j 

The  mean  compliance  of  the  mantle  ( J  =  l//x) 

j 

When  there  is  no  risk  of  confusion,  J  also  denotes 
the  Jacobian 

p 

The  mean  density  of  the  planet 

g 

The  surface  gravity  on  the  planet 

u 

Tidal  displacement  in  the  planet 

V 

The  velocity  of  tidal  displacement  in  the  planet 

p 

The  power  exerted  on  the  planet  by  the  tidal  stresses 

( p ) 

The  time- averaged  power  (averaged  over  one  or 
several  cycles  of  flexure) 

2.3.  Difficulties 


We  shall  derive  the  heat-production  formulae  for  two  dif¬ 
ferent  settings — with  a  fixed  pericenter  to  and  with  o>  moving 
uniformly. 


2.4.  Lagging 

For  a  static  tide,  the  incremental  tidal  potential  of  the 
perturbed  body  mimics  the  perturbation  (Equation  (2)),  except 
that  each  term  Wi  is  now  equipped  with  a  mitigating  multiplier 
^7  (R/r)l+l,  where  k /  is  an  /-degree  Love  number.  With  the  star 
located  in  r  *,  the  additional  potential  in  a  point  r  will  read  as 


U(r,  r  *)  =  £t/,(r) 


1=2 


-EM  t 


1=2 


R 


l+\ 


W,(R.  /•*). 


(10) 


At  this  point,  a  word  of  warning  is  necessary.  Deriving 
Equation  (5),  we  differentiated  the  expression  (3),  which  gave 
us  the  terms  with  to  and  LI.  Including  these  terms  in  Equation  (5) 
acknowledges  the  fact  that  the  perturber’s  trajectory  is  disturbed, 
not  Keplerian.  The  disturbance  may  come  solely  from  tides,  as 
in  Kaula  (1964,  Equation  (38)),  or  from  both  tides  and  other 
sources.  One  way  or  another,  the  perturber’s  mean  anomaly  M 
is  no  longer  equal  to  it,  but  is  now  given  by  M  =  Mq  + 
f*  n(t)dt ,  where  n{t)  =  2^44MDMD)a~Htj .  Accordingly, 
the  expression  for  the  modes  becomes: 

lOlmpq  —  illmpq  171  0  =  (/  2 p)  (O)  +  MIq) 

+  q  Mo  +  (/  —  2 p  +  q)n  +  m  (LI  —  0).  (7) 

It  is,  of  course,  tempting  to  assume  that  Mo  n,  thus  accepting 

the  approximation 


For  time-dependent  tides,  this  expression  acquires  an  extra 
amendment:  the  reaction  must  lag,  compared  to  the  action. 
Naively,  this  would  imply  taking  each  W/  at  an  earlier  instant  of 
time.  However,  in  reality  lagging  depends  on  frequency,  so  each 
W/  must  be  first  expanded  into  a  Fourier  series  over  tidal  modes, 
whereafter  each  term  of  the  series  should  be  delayed  separately. 
The  magnitude  of  the  tidal  reaction  is  frequency  dependent 
too,  so  each  term  of  the  Fourier  series  will  be  multiplied  by 
a  dynamical  Love  number  of  its  own.  Symbolically,  this  may  be 
written  in  a  manner  similar  to  the  static  expression: 


U(r,  r*)=  ^2  U/(r ,  r*) 


1=2 


=  £ 


1=2 


R 


i+ 1 


k,  W,(R.  r  *). 


(11) 


M  =  Mo  +  n  ^  n,  (8) 

as  Kaula  (1964)  did  in  his  Equations  (46)  and  (47).  Within  his 
theory,  however,  this  approximation  could  not  be  used.4  This 
is  explained  in  Appendix  B  where  we  consider  two  examples. 
One  is  the  case  where  perturbation  of  an  orbit  of  a  moon  is 
mainly  due  to  the  tides  the  moon  creates  in  the  planet.  In  that 
situation,  to  and  Mo  are  of  the  same  order  but  of  opposite 
signs,  so  they  largely  compensate  one  another.  This  suggests  a 
simultaneous  neglect  of  both  rates.  The  second  example  is  when 
the  dominant  perturbation  of  the  orbit  comes  from  the  oblateness 
of  the  primary.  In  this  case  to  and  Mo  are  of  the  same  order  and 
the  same  sign — so  keeping  one  of  these  terms  requires  keeping 
the  other. 

Whether  one  or  both  of  these  rates  should  be  included  depends 
on  a  particular  setting,  and  each  practical  case  must  be  examined 
separately.  In  general,  both  rates  should  be  kept. 

While  keeping  to  complicates  the  formalism,  the  emergence 
of  Mo  complicates  the  treatment  even  further.  To  sidestep  this 
issue,  we  shall  define  the  mean  motion  via 


The  hat  above  ki  means  that  this  is  not  a  multiplier  but  a  linear 
operator  that  mitigates  and  delays  each  Fourier  mode  of  W/ 
differently: 


U(r, 


a 


OO 


£ 


/ 

£ 

m= 0 


(/  —  m)\ 
(/  +  m)\ 


X  (2  0)  ^  ^  ^'impH') 

P- 0 


oo 

X  ^  ^  Glpq{c)kl{(L)lmpq ) 

q=— oo 

X  ( Vlmpq  -  m  (X  +  0)  -  ei ), 


cos 

sin 


l — III  even 


(12) 


where  the  Love  numbers  k/(toimpq )  and  the  phase  lags  ei(a>impq) 
are  functions  of  the  Fourier  modes.  The  lags  emerge  as  the 
products 


C  /( (lllmpq)  —  tOlmpqtSti(,tO\mpq ),  (13a) 


n  =  M.  (9) 

This,  the  so-called  anomalistic  mean  motion  differs  from 
JG(M  +  M*)a~\t). 


4  In  his  books,  Kaula  (1966,  1968)  corrected  this  oversight.  There,  he  kept 
the  notation  n  for  the  mean  motion  defined  as  in  the  Kepler  law,  and  never 
confused  it  with  M  . 


where  A ti(coimpq)  is  the  time  delay  at  the  mode  toimpq.  In  reality, 
the  time  delays  are  functions  not  of  the  Fourier  modes  (which 
can  assume  either  sign),  but  of  the  actual  physical  forcing 
frequencies  ximpq  —  which  are  positive  definite.  Thus 

it  is  more  accurate  to  write  the  delays  not  as  A ti(toimpq)  but  as 
A ti(ximpq).  Accordingly,  the  phase  lags  become 

tlitOlmpq )  =  tOlnipq  Atfix  Impq  )  -  (13b) 
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For  the  reasons  of  causality,  the  time  delays  are  positive  definite, 
so  the  sign  of  the  phase  lag  always  coincides  with  that  of  the 
corresponding  Fourier  mode.  Thus  we  finally  have: 

e/ i^lmpq)  —  (*)lmpq  Sgll  (o)ln!pq)  Nt\fXlmpq) 

—  Xlmpq  SgU  i^lmpq)  ^tfXlmpqf  (13c) 

where  xtmpq  =  \u>impq\  are  the  positive  definite  forcing 
frequencies. 

The  dynamical  Love  number  ki(coimpq)  and  the  phase  lag 
ei(ojimpq)  are  the  absolute  value  and  the  negative  phase  of  the 
complex  Love  number  ki(a>impq )  whose  functional  dependence 
upon  the  Fourier  mode  is  solely  determined  by  /,  provided  the 
body  is  spherical.' 

2.5.  Physics  behind  the  Love  Numbers  and  Phase  Lags 

As  we  saw  above,  to  obtain  the  decomposition  (12)  from  the 
Fourier  series  (2),  each  Impq  term  of  the  latter  had  to  be  endowed 
with  its  own  mitigating  factor  k /  =  k[(a>impq )  and  phase  lag 
6/  =  €i(ojimpq).  In  the  past,  some  authors  enquired  whether  this 
mitigate-and-lag  method  is  general  enough  to  describe  tides.  It 
is,  as  long  as  the  tides  are  linear.  This  is  explained  in  Appendix  C. 

The  expression  (12)  for  the  additional  tidal  potential  contains 
both  sines  and  cosines  of  the  phase  lags,  and  so  does  the  ensuing 
expression  for  the  surface  elevation.  However,  the  resulting 
expression  for  the  tidal  dissipation  rate  turns  out  to  contain  only 
the  combination  kfco)  sin  €/(&>)  which  is  the  negative  imaginary 
part  of  the  complex  Love  number: 

ki(a>)  sine/(tt>)=  |£/(a>)|  sine/(<a)=  —  lm[ki(a>)], 
where  co  =  u>hnpq .  (14) 

This  quantity  is  often  denoted  as  ki/Q ,  although  it  would  be 
more  reasonable  to  employ  the  notation  kfi  Q\,  with  the  tidal 
quality  factors  defined  through  l/Qi  =  sin  f/|. 

A  dynamical  Love  number  ki(a>impq)  is  an  even  function  of 
the  tidal  mode  a>impq,  while  a  phase  lag  ei(coimpq)  is  odd,  as 
can  be  observed  from  Equation  (13b).  Thus  the  expression  for 
the  product  k\  sine/  as  a  function  of  the  physical  frequency 

X  ~  Xlmpq  —  t^lmpq  IS. 

k,{co)  sin €i(co)  =  kfx)  sine/(x)  Sgn(«),  (15) 

where  €/(/)  is  non-negative,  because  the  physical  frequency 
is  x- 

The  frequency  dependence  of  k//Qi  =  ki(x)  sine/(x)  is 
defined  by  two  major  physical  circumstances:  self-gravitation 
of  the  planet  and  the  rheology  of  its  mantle.  A  rheological  law 
is  expressed  by  a  constitutive  equation,  that  is,  by  an  equation 
interconnecting  the  strain  and  the  stress.  A  particular  form  of  this 
equation  is  determined  by  the  friction  mechanisms  present  in  the 
considered  medium.  A  realistic  rheological  law  should  contain 
contributions  from  elasticity,  viscosity,  and  inelastic  processes 
(mainly,  dislocation  unjamming).  Self-gravitation  suppresses 
the  tidal  bulge.  At  low  frequencies  this  effectively  adds  to  the 
mantle’s  rigidity,  whereas  at  higher  frequencies  the  interplay 
of  rheology  and  gravity  is  more  complex  (Efroimsky  2012b, 
Figure  2). 


5  For  oblate  celestial  bodies,  the  functional  form  of  the  complex  ki(u>impq)  is 
also  determined  by  the  order  m.  In  that  situation,  the  right  notation  for  the 
complex  Love  number  is:  kim(a>impq).  Its  absolute  value  and  negative  phase 
will  then  be  denoted  with  kim(coimpq)  and  tim(u>impq). 


The  calculation  of  the  frequency  dependence  kfx)  sin  <?/(x) 
for  a  homogeneous  body  of  a  known  size,  mass,  and  rheology 
is  presented  in  detail  in  Efroimsky  (2012a,  2012b).  See  also  the 
Appendix  to  Paper  II. 

While  quadrupole  (/  =  2)  terms  are  sufficient  in  most 
problems,  exceptions  are  known.  For  the  orbital  evolution  of 
Phobos,  the  /  =  3  and,  possibly,  even  1  =  4  terms  of  the 
Martian  tidal  potential  may  be  of  relevance  (Bills  et  al.  2005). 
Studying  close  binary  asteroids,  Taylor  &  Margot  (2010)  took 
into  account  the  Love  numbers  up  to  /  =  6. 

The  question  of  how  rapidly  /  >  2  terms  fall  off  with  the 
increase  of  the  degree  /  is  also  interesting.  Most  authors  only 
rely  on  the  geometric  factor  ( R/a)2l+l  to  answer  this  question. 
As  was  explained  in  Efroimsky  (2012b),  the  /-dependence  of 
ki(a>impq)  sin  e(a>impq),  too,  comes  into  play  and  changes  the 
result  considerably. 

3.  THE  EULERIAN  AND  LAGRANGIAN  DESCRIPTIONS 

What  we  hope  ever  to  do  with  ease, 
we  must  learn  first  to  do  with  diligence. 

Samuel  Johnson 

3.1.  Notations  and  Definitions 

To  compare  the  varying  shape  of  a  deformable  body  against 
some  benchmark  configuration,  we  use  X  to  denote  the  initial 
position  occupied  by  a  particle  at  t  =  0.  At  another  time  f,  the 
particle  finds  itself  in  a  new  place 

*  =  f{X,  t), 

where  the  function  f(X,t)  is  a  trajectory,  that  is,  a  solution  to 
the  equation  of  motion,  with  the  initial  condition  x  =  X  set  at 
t  =  0. 

The  current  values  of  all  physical  and  kinematic  properties  of 
the  medium  can  be  expressed  as  functions  of  the  instantaneous 
coordinates  x  of  a  point  where  these  properties  are  being 
measured  at  the  present  moment  t.  When  referred  to  the  present 
time  and  position,  such  properties  are  named  Eulerian  and 
are  equipped  with  a  subscript  E\  for  example:  qE(x.  t).  The 
Eulerian  description  is  fit  to  answer  the  question  “where,”  and 
therefore  is  convenient  in  fluid  dynamics  where  the  displacement 
x  —  X  can  become  arbitrarily  large  and  the  initial  position  X  is 
soon  forgotten. 

While  x  denotes  a  place  in  space,  the  initial  condition  X  acts 
as  the  “number”  of  a  particle  presently  residing  at  the  place  x. 
Although  located  at  jc,  the  particle  originally  came  from  X  and 
will  carry  the  label  X  forever. 

Knowing  the  trajectories  of  all  particles,  we  can  express  the 
properties  as  functions  of  the  time  t  and  the  initial  conditions 
X.  To  that  end,  we  employ  the  change  of  variables  x  = 
f(X,  t ).  Expressed  through  the  initial  conditions,  a  property 
q  will  be  termed  as  Lagrangian  and  equipped  with  the  subscript 
L: 

qL(X,  t)  =  qE(x ,  t)  (16a) 

or,  in  more  detail: 

qL{X,t)  =  qE{f{X,t),t).  (16b) 

So  qE  has  the  same  value  as  q E ,  but  has  a  different  functional 
form,  as  it  is  now  understood  as  a  function  of  the  initial 
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conditions  (the  particles’  “numbers”)  X,  and  not  of  the  present¬ 
time  coordinates  x.  Relating  the  quantities  to  the  initial  positions 
X ,  the  Lagrangian  description  tells  us  “which  particle”,  and  is 
thus  practical  in  description  of  deformable  solids. 

In  anticipation  of  perturbative  treatment,  we  regard  the 
trajectory  x  =  f(X,  t)  as  fiducial  and  equip  the  appropriate 
functional  dependencies  with  a  superscript  0: 

q°L(X,  t )  =  q°E(x,  t)  (17a) 

which  is: 

q°L(X,t)  =  q°E(f(X,t),t).  (17b) 

3.2.  Perturbative  Approach 

Under  disturbance,  two  changes  will  take  place  in  a  point  r 
at  a  time  t: 

1.  Properties  will  now  assume  different  values  in  this  point 
at  this  time.  So  we  substitute  the  unperturbed  Eulerian 
dependencies  qE(r,  t )  with 

qE(r,t)  =  qE(r,t)  +  q'E(r,t).  (18) 

This  equality,  in  fact,  serves  as  a  definition  of  the  variation 
q'E{r,  t ):  the  variation  is  a  change  in  the  functional  depen¬ 
dence  of  a  physical  property  upon  the  present  position  r. 

2.  A  different  particle  will  now  appear  in  the  point  r  at  the 
time  t.  It  will  not  be  the  same  particle  as  the  one  expected 
there  at  the  time  t  in  the  absence  of  perturbation. 

Accordingly,  a  particle,  which  starts  in  X  at  t  =  0,  will 
show  up  at  the  time  t,  not  in  the  point  x  —  f(X,  t )  but  in 
some  other  location  displaced  by  u: 

r  =  x+u=  f(X,t)  +  u(X,t).  (19) 

Both  of  these  changes,  1  and  2,  will  affect  the  Lagrangian 
dependencies  of  the  properties  upon  the  initial  conditions,  so  the 
dependency  of  each  property  will  acquire  a  variation  q'L(X .  t): 

qL(X,t)  =  q0L(X,t)  +  q'L(X,t).  (20) 

In  Appendix  D,  we  provide  a  self-sufficient  introduction  into  the 
perturbative  treatment  of  a  deformable  body,  both  in  the  Eulerian 
and  Lagrangian  languages.  There  we  derive  a  relation  between 
the  perturbations  of  the  Lagrangian  and  Eulerian  quantities: 

q'L(X,  t )  =  q'E(f{X,  t),  t)  +  u(X,  t)VxqE  +  0( u2),  (21) 

with  the  gradient  in  the  second  term  acting  on  the  unperturbed 
history:*1 

VxqE  =  VxqE(x,  t ),  where  x  =  f(X ,  t ).  (22) 

In  the  formula  (21),  the  first  term  on  the  right-hand  side,  q'E(x,  t), 
accounts  for  the  change  of  the  final  spatial  distribution  of 
properties.  The  other  two  terms  show  up  because  perturbation 
alters  the  mapping  from  X  to  the  present  position. 


6  To  derive  Equation  (21),  we  expanded  qE(r,  t)  =  qE{x  +  u,  /)  into  the 
Taylor  series  near  the  unperturbed  qE(x,  t). 


3.3.  Summary  of  Linearized  Formulae  for  the  Density  of  a 
Periodically  Deformed  Solid 

We  need  several  formulae  for  density  perturbations,  which 
are  obtained  in  Appendix  D. 

In  the  Eulerian  description: 


pE(r,  t )  =  pE(r)  +  p'E(r,  t ), 

(23) 

E*J  " 

+ 

<1 

■'S 

ft]  o 

II 

O 

(24) 

Formula  (23)  renders  the  interrelation  between  the  functions 
of  the  same  variable.  The  unperturbed  density  pE  appears  here 
as  a  function  of  the  perturbed  present  positions  r,  not  of  the 
unperturbed  reference  positions  x.  This  can  be  traced  through 
the  derivation  (D21)-(D24).  There,  the  unperturbed  density 
initially  shows  up  as  a  function  of  x  —  r  —  u.  It  then  ends 
up  as  a  function  of  r,  after  the  Taylor  expansion  around  r  over 
powers  of  u  is  performed. 

Accordingly,  the  symbol  Vr  denotes  differentiation  with 
respect  to  the  perturbed  position  r  upon  which  pE  is  set 
to  depend  in  the  above  equations.  Also  remember  that  in 
u(x,  t)  =  u(r,  t)  +  (9(u2)  we  can  neglect  0(u2),  in  the  linear 
approximation.  Thus  the  Lagrangian  and  Eulerian  values  of 
the  displacement  coincide  in  the  first  order.  Specifically,  in 
Equation  (24),  our  u  can  be  treated  as  a  function  of  r.  So  all 
entities  in  that  equation  are  functions  of  the  same  variable,  the 
perturbed  location. 

In  the  Lagrangian  description: 

pL{X,t)  =  p°E{X)  +  p'L{X,t),  (25) 

p'L  +  Pe  Vx  •  u  =  0.  (26) 

Recall  that  this  is  an  interrelation  between  functions  of  the 
same  variable.  This  time,  it  is  the  initial  position  X.  Had  we 
altered  the  notation  from  T  to  r,  nothing  would  have  changed 
(except  that  we  would  write  V,.  instead  of  Vy ) — it  would  still  be 
the  same  relation  between  three  functions  of  the  same  argument. 

Relation  between  the  increments  p'r  and  p'F  \ 

This  relation  originates  from  the  general  formula  (21).  In  our 
case,  the  reference  trajectory  x  —  ,f(X.  t )  stays  identical  to  the 
initial  position  X,  so  we  obtain: 

p'L(x,  t)  =  p'E(x,  t)  +  u(x,  t)  •  V,  p°(x,  t).  (27) 

Once  again,  we  are  dealing  with  a  relation  between  several 
functions  taken  all  at  one  and  the  same  point.  Here  the  point 
is  denoted  with  x.  Had  we  denoted  it  with  r,  the  only  change 
would  be  a  switch  from  Vv  to  Vr,  no  matter  what  meaning  we 
instill  into  these  x  and  r. 

For  an  initially  homogeneous  body,  V  pE  =  0,  so  the 
forms  (24)  and  (26)  of  the  linearized  conservation  law  coincide 
and  can  both  be  conveniently  written  as 

p°Vf. -v  +  ^  =0,  (28) 

ot 

where  p  0  =  pE  and  the  velocity  is 
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3.4.  Potentials  and  Their  Increments 

In  each  point,  the  density  p  and  potential  V  comprise  a  mean 
value  and  a  perturbation: 


the  integration  being  performed  over  an  instantaneous,  deformed 
volume.  Together  with 

p£v  •  V,.  V'E  =  V,.  •  {pE  yV’E)  -  V'E  V,.  •  (pE  v),  (37) 


density:  p  =  p°  +  p',  (30a) 

potential:  V  =  V  0  +  V'  =  V°  +  (W  +  U),  (30b) 

where  V  0  is  the  constant-in-time  spherically  symmetrical  po¬ 
tential  of  an  undeformed  body,  whereas  V'  denotes  the  po¬ 
tential’s  perturbation.  The  perturbation  consists  of  the  external 
tide-raising  potential  W  and  the  resulting  additional  potential  U 
of  the  perturbed  body: 


V'=W  +  U.  (31) 

The  potentials  and  densities  will  be  endowed  with  a  subscript 
“L  ”  or  “E"  pointing  at  the  Lagrangian  or  Eulerian  descriptions, 
accordingly.  Owing  to  the  general  expression  (21),  we  have: 

VE(r,  t)  =  VL'(x,  0  -  u  ■  Vx  VE,  (32) 

the  same  being  valid  for  p,  see  Equation  (27).  For  unperturbed 
properties,  however,  subscripts  may  be  dropped  without  causing 
any  confusion: 


Va=V°,  p°  =  p°.  (33) 

3.5.  The  Poisson  Equation  in  the  Eulerian  Description 

In  both  the  perturbed  and  unperturbed  settings,  the  density 
and  potential  are  always  linked  through  the  Poisson  equation: 

V,.2  VE  =  -4?r  GPe,  (34a) 

V,.2  VE  =  —An  G  pE,  (34b) 


the  mass-conservation  law 

V,-  -(p£v)+  ^  =  0  (38) 

at 

simplifies  the  expression  under  the  integral  to  the  following 
form: 

p£v  ■  Vr  V'E  =  Vr  ■  (pE  v  V'E)  +  VE'^f-  .  (39a) 

dt 

Further  employment  of  the  Poisson  equation  in  the  Eulerian 
form,  (Equation  (35)),  gives  us 

1  3 

pEy  ■  Vr  V’,  =  V,.  ■  (pE  v  V'E)  -  —  V'E  -  V2  VE. 

(39b) 

So  the  power  becomes 

P=  J  V,.  ■  (pE  v  V'E)  d\  +  J  V'E^d^r  (40a) 

=  [  pEVE'\-dS 1 
Jr 

(40b) 

where  d S'  =  nf  dY! ,  with  n'  and  dV  being  a  unit  normal  to  the 
deformed  surface  and  an  element  of  area  on  that  surface,  both 
taken  at  the  time  t.  Correct  to  the  first  order  in  the  displacement 
u,  these  are  related  to  their  unperturbed  analogues  via 


n'  =  (1  —  Vz  <g>  u)n  and  dl!  =  (1  +  Vz  •  u)dY, 

(41) 


where  the  surface  gradient  is  defined  as 


while  the  perturbing  potential  W  obeys  the  Laplace  equation 
outside  the  perturber: 

V2  WE  =  0.  (34c) 

Subtraction  of  Equation  (34b)  from  Equation  (34a)  results  in 
a  Poisson  equation  for  the  density  perturbation: 

Vr2  V'  =  -4nGp'E.  (35) 

The  Poisson  equation  in  the  Lagrangian  description  is  presented 
in  Appendix  D. 

4.  THE  POWER  PRODUCED  BY  THE  TIDAL  FORCE 

4.1.  In  the  Eulerian  Description 

The  power  P  exerted  on  the  perturbed  body  is  an  integral, 
over  its  volume,  of  the  rate  of  working  by  tidal  forces  on 
displacements.  In  the  Eulerian  language,  the  power  reads  as 


P  =  I  PeV  ■  Vr  V'E  d\,  (36) 


vz  =  Vx  -  n9fl,  (42) 

so  V1  <g>  u  is  a  three-dimensional,  second-rank  tensor  (Dahlen 
&  Tromp  1998).  Altogether, 

dS'  =  n'  dY!  =  (1  +  Vz  •  u)  MY  -  (Vz  <g>  u)  MY 

=  (1+  Vz-uWS-(Vs®u)JS,  (43a) 

with  d S  =  MY  pertaining  to  the  unperturbed  surface.  In  a 

shorter  form,  the  above  reads  as 

d&  =  IdS,  (43b) 

where  the  three-dimensional,  second-rank  tensor 

Jl  =  (l+  Vz-u)I-  Vs  ®u  (44) 

is,  loosely  speaking,  playing  the  role  of  a  Jacobian  for  elements 
of  area.  This  is  fully  analogous  to  the  formula 

cPr  =  J  d^x  =  (1  +  V*  •  u)d3x  =  [  1  +  V,.  •  u  +  0(u2)]d3x 

(45) 

linking  the  deformed  volume  d3r  to  the  undeformed  volume 
d3x  (see  Appendix  DAL). 
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4.2.  In  the  Lagrangian  Description 
Applied  to  the  density,  the  general  formula  (16)  renders: 


pL(X,  t)  =  pE(r ,  t).  (46) 

This,  together  with  the  formula  (32)  for  the  potential  pertur¬ 
bation,  enables  us  to  express  the  power  in  the  Lagrangian 
description: 

P  —  J  Pl  v  •  V,  VL'  d3x  -  J  pL  v  ■  V*  (u  •  V.t  V°)  d3x, 

(47) 

the  integral  now  being  taken  over  the  undeformed  body.  Be 
mindful  that  d3r  Vr  =  d3x  Vt ,  so  no  Jacobian  shows  up  on  the 
right-hand  side. 

The  velocity  and  displacement  being  in  quadrature,  the 
second  term  should  be  dropped  after  time  averaging  (denoted 
with  angular  brackets): 

(P)  —  J  Pl  v  •  V,  Vl  d3x,  (48a) 

For  a  periodically  deformed  solid,  we  set  the  equilibrium  state 
to  play  the  role  of  the  unperturbed  configuration,  for  which 
reason !  pE( X,  t)J  —  pE(x).  Insertion  of  this  equality  into  the 
expression  (48a)  gives  us: 

(P)  =  J  p°v  •  VXVL'  J~1 *d3x.  (48b) 

The  dot-product  can  be  easily  rearranged  via  the  formulae 
analogous  to  Equations  (37) — (39).  Due  to 

P °v  •  V*  VL'  =  V*  •  (p  0  v  VL')  -  VL'  V,  •  (v  p  °)  (49) 

and 

n  9p 

P°vy-v+/=°,  (50) 

ot 

the  expression  under  the  integral  becomes 

p°y  ■  V,  V'L  =  V,  •  (p°  v  V'L)  +  VL'^,  (51) 

provided  we  set  Vv p°  —  0,  that  is,  provided  we  assume  that 
the  unperturbed  body  is  homogeneous.8  Then  the  time-averaged 
power,  for  an  initially  homogeneous  body,  acquires  the  form  of 

( P )  =  J  V*  •  (p°  v  VL')  d3  x  +  J  VL'^d3x,  (52) 

where  we  approximated  the  Jacobian  with  unity,  thus  neglecting 
higher-order  terms. 

1  The  mass  is  conserved  along  both  trajectories,  perturbed  and  unperturbed. 

So  both  pE(r,  t)d3r  and  pE(x,  t)d^x  must  be  equal  to  the  initial  mass 

pE(X)  d3 X ,  and  therefore  to  one  another:  pE(r,  t)d^r  =  paE(x,  t)d3x. 

Thence,  pE(r,  t)J  =  pE(X),  where  J  =  d?r /d?x.  In  combination  with 

Equation  (46),  this  yields:  pE(X ,  t)J  =  pE(x.  t). 

When  the  unperturbed  configuration  is  the  equilibrium  state,  x  =  f(X.  t) 

coincides  with  X  at  all  times.  So  pE(x ,  f)  bears  no  dependence  on  time,  and 
the  above  equality  becomes  simply  pE (X,  t)J  =  pE(x).  See  Appendix  D.4.1 
for  a  detailed  discussion. 

8  No  such  assumption  was  required  to  obtain  the  Eulerian  analogue 
(Equation  (38))  of  the  Lagrangian  formula  (50). 


5.  TIDAL  DISSIPATION  RATE 
IN  A  HOMOGENEOUS  SPHERE 

Although  the  Eulerian  and  Lagrangian  pictures  are  equivalent, 
the  boundary  conditions  look  simpler  in  the  Eulerian  version.  On 
the  other  hand,  for  periodic  deformations,  practical  calculations 
are  easier  carried  out  in  the  Lagrangian  description,  as  it 
implies  integrations  over  the  unperturbed  volume  and  surface 
corresponding  to  the  equilibrium  shape.  It  is,  unfortunately, 
not  unusual  for  the  authors  to  refrain  from  pointing  out  which 
description  is  employed,  leaving  this  to  the  discernment  of  the 
readers.  The  easiest  way  to  trace  an  author’s  choice  is  to  look  at 
the  way  they  write  the  expression  for  the  power  and  the  Poisson 
equation. 

The  often-cited  authors  Zschau  (1978)  and  Platzman  (1984) 
started  in  the  Eulerian  language  and  then  switched  to  the 
Lagrangian  description.  This  can  be  seen  from  the  fact  that  the 
time-average  power  was  eventually  written  by  both  of  them  as  an 
integral  over  the  undefonned  body.  Both  works  contained  some 
mathematical  omissions,  which,  fortunately,  did  not  influence 
the  final  form  of  the  integral. 

Below  we  present  these  authors’  method  in  a  more  mathe¬ 
matically  complete  manner.  While  our  expression  for  the  power, 
written  as  an  integral  over  the  unperturbed  surface,  will  coincide 
with  the  integrals  derived  by  the  said  authors,  our  final  result 
(the  power  written  as  a  spectral  sum  over  the  Fourier  modes) 
will  differ.  In  one  important  detail,  our  result  also  differs  from 
that  by  Peale  &  Cassen  (1978). 


5.7.  A  Mixed,  Eulerian— Lagrangian  Treatment 

Similar  to  Zschau  (1978,  Equation  (2)),  we  begin  with  the 
formula  (36)  for  the  power  in  the  Eulerian  variables.  The  next 
natural  step  is  Equation  (40),  whereafter  integration  by  parts 
renders: 


F=  jPEV^.dS -A- 


L  [ 

4itG  J 


3  avrv£' 


d'r 


dt 


vr  V 


9Vr  Ve'  \ 

dt  J 


(53a) 


f  If  ,dXrVE' 

=  J  pLVE'y-QdS)-  —  J  VE'— r^~  -QdS) 

13/* 

+  Ve'  ■  Vr  Ve'.  (53b) 

sjtG  dt  J 

En  route  from  the  former  expression  to  the  latter,  we  switch 
from  dS‘  and  d3r  to  J  <7S  and  J  d3x,  respectively.  Thereby 
we  switch  from  integration  over  a  deformed  body  to  that  over 
the  undeformed  one.  So  pE  becomes  pE,  see  Equation  (46). 
A  similar  switch  from  VE  to  VE  can  be  performed  using 
Equation  (32),  but  we  prefer  to  stick  to  VE  for  some  time, 
as  it  will  be  easier  to  impose  the  boundary  conditions  on  the 
Eulerian  potential. 

In  a  leading-order  calculation,  both  the  Jacobian  and  its 
tensorial  analogue  may  be  set  unity:  J  ^  I  and  J  ~  1,  as 
evident  from  the  formulae  (44)  and  (45).  In  the  same  order, 
we  can  substitute  V,.  with  Vr .  In  addition,  as  was  explained  in 
Footnote  7,  we  can  substitute  pE  —  p°/J  with  p°,  and  treat  the 
latter  as  time-independent.  Thus,  the  time  average  of  the  power 
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( P ) 


/ 


(P°Ve'v)  •  d S- 


3  VXVE' 
3  t 


d  S 
(54a) 


(54b) 

with  the  volume  integral  dropped. The  potential  VE  in  the 
above  developments  was  the  interior  potential,  so  the  above 
formula  should,  rigorously  speaking,  have  been  written  as 

/p\  = _ ! _  f  /  y  '(interior)  /vy  y  /(interior) 

'  4nG  J  \  dtK 

—  4jt  G  p°u) )  •  dS.  (54c) 


due  to  luck  than  to  accuracy.  In  the  subsequent  derivation,  the 
author’s  formulae  (7)  and  (10)  are  incorrect  because  the  fact  that 
the  Fourier  modes  in  the  Darwin-Kaula  theory  can  be  of  either 
sign  is  neglected.  We  address  this  point  at  the  end  of  Section  5.4. 

5.3.  Employment  of  the  Boundary  Conditions 

The  Eulerian  boundary  conditions  mimic  those  from  electro¬ 
statics  (see  Appendix  E): 


(interior)  (exterior) 

VE'  =  VE' 


(55) 


and 


(exterior) 


—  VE'  -  4tt  Gp  °u 

3n 


—  VE'-47tGp°u 

3n 


(interior) 

.  (56) 


Insertion  thereof  into  Equation  (54c)  makes  the  power  look 


The  expression  (54c)  is  somewhat  formal.  On  the  one  hand, 
it  contains  integration  over  an  undeformed  surface,  an  operation 
appropriate  to  the  Lagrangian  description.  On  the  other  hand, 
the  quantity  under  the  integral  is  Eulerian,  that  is,  a  function  of 
the  perturbed  positions.  Thus,  to  employ  the  expression  (54c) 
in  practical  calculations,  one  would  first  have  to  express  the 
integrated  average  product  (V£,<lntenor,(3/3f)  (Vt  Ve'  (mtenol)  — 
4  it  G  f) 0  li  J )  as  a  function  of  the  unperturbed  positions,  that  is, 
of  the  coordinates  on  the  undeformed  surface.  Simply  speaking, 
one  would  have  to  switch  from  a  Eulerian  function  under  the 
integral  to  a  Lagrangian  function,  using  the  formula  (32).  The 
reason  for  our  procrastination  with  this  step  is  the  convenience 
of  the  Eulerian  description  for  imposing  boundary  conditions. 

5.2.  Comparing  the  Intermediate  Result  (Equation  (54c))  with 
Analogous  Formulae  from  Zschau  (1978)  and  Platzman  (1984) 

Our  expression  (54c)  is  equivalent  to  formula  (12)  in  Zschau 
(1978).  The  sole  difference  is  how  we  justify  the  substitution  of 
the  Lagrangian  density  pE  with  the  unperturbed  p°.  Whereas 
we  approximated  the  Jacobian  with  1  +  0(|u|),  Zschau  (1978, 
Equation  (10))  employed  a  clever  trick  that  did  not  rely  on  the 
smallness  of  disturbance.  In  our  notation,  the  trick  looks  like 
this:  if  in  the  first  term  of  our  expression  (54a)  we  also  keep  the 
first-order  perturbation  p’L  of  the  density,  the  time  average  of 
the  product  p'L  v  VE'  will  always  be  zero,  provided  all  three 
oscillate  at  the  same  frequency.  While  elegant,  Zschau’s  argu¬ 
ment  works  only  for  a  perturbation  at  one  frequency,  not  for  a 
spectrum  of  frequencies. 

The  treatment  by  Platzman  (1984)  contains  more  inaccura¬ 
cies.  The  author’s  formula  (2)  looks  like  our  Equation  (48b), 
with  the  actual  density  substituted  from  the  beginning  by  its 
unperturbed  value  p°.  Such  a  start  indicates  the  use  of  the  La¬ 
grangian  description.  This  however  comes  into  contradiction 
with  the  way  the  author  writes  the  conservation  law.  Platzman’s 
form  of  that  law  is  equivalent  to  our  Equation  (24) — it  is  writ¬ 
ten  in  the  Eulerian  language.  The  following  Poisson  equation 
is  also  Eulerian.  That  the  author  eventually  arrives  at  the  right 
integral  expression  (Equation  (5)  in  Platzman  (1984))  is  more 


9  As  previously  agreed,  in  our  approximation  the  Jacobian  is  set  to  unity.  The 
potential  variation  Ve '  is  a  sum  of  sinusoidal  harmonics,  and  so  is  its  gradient 
Vx  Ve' ■  After  time  averaging  of  the  expression  (53b),  the  cross  terms  in  the 
product  Vx  Ve'  ■  VA  V p '  will  vanish,  while  the  products  of  harmonics  of  the 
same  frequency  will  render  constants. 


(P)  =  ~ 


4jtG 


T/  r  (exterior) 

Ve 


3 

3 7 


17  /  (exterior) 
V  x  V  E 


d  S.  (57) 


It  is  now  time  to  write  the  expression  under  the  integral  (57) 
as  a  function  of  the  coordinates  on  the  unperturbed  surface,  the 
one  over  which  we  integrate.  The  formula  (32)  prescribes  us  to 
substitute  VE  with  Vjf  —  u  •  VA  Vo.  As  u  is  zero  outside  the  body, 
we  get10 


(P) 


=  ~—f 

4nG  J 


VL 


(exterior)  3 


—  V  V  ,fexterior) 
3 1  x  L 


■  d S.  (58) 


To  analyze  the  behavior  of  V'  outside  the  perturbed  body, 
recall  that  its  two  components,  U  and  W,  scale  differently  with 
the  planetocentric  radius.  As  can  be  seen  from  Equation  (1), 
the  degree-/  Legendre  component  of  the  perturbing  potential 
changes  as11  IV/  oc  r1 .  According  to  Equation  (II),  the  degree- 
/  component  of  the  tidal  potential  obeys  17/  oc  r  (,+  l).  All  in  all, 
the  /-degree  part  of  the  exterior  V'  assumes  the  form  of 


VL' 


(exterior) 


-(7+1) 


Ui(R) 


1=2  L 


(59) 


while  the  normal  part  of  its  gradient  on  the  free  surface  is 


3 

dr 


(exterior) 

Vl', 


=  /r1 


(iw,  -  (/  +  \)U,  ] . 

1=2 


(60) 


10  For  the  first  multiplier  under  the  integral  (57),  we  simply  substitute 

Ve'  (extenor)  wjth  yL  f  (exterior) ^  omjttjng  term  [_u  .  Vx  y0]  (extenor)  because  U 
is  zero  outside  the  body.  The  case  of  the  second  multiplier, 

(d/dt)Vx  js  jess  obvious.  Employment  of  the  formula  (32)  gives 

(3/3 t)[Vx  Vl  —  V*  •  (u  V  °)]  (extenor).  The  vanishing  of  u  on  the  exterior  side 
of  the  boundary  does  not  imply  the  vanishing  of  its  gradient  there.  On  the 
contrary,  V*  •  (u  V  °)  performs  a  finite  step,  but  so  also  does  the  gradient  of 
V^',  so  that  altogether  the  gradient  Ve'  remains  continuous.  To  sidestep  these 
intricacies,  we  can  expand  the  volume  of  integration  slightly  outward  from  the 
actual  volume  of  the  planet  (Platzman  1984,  p.  74). 

11  Do  not  be  misled  by  the  planetocentric  distance  in  Equation  (1)  being 
denoted  with  R.  There  we  needed  the  value  of  W  on  the  surface,  whereas  here 
we  need  to  know  W  at  an  arbitrary  planetocentric  distance. 
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Plugging  it  into  Equation  (58),  and  benefiting  from  the  orthog¬ 
onality  of  surface  harmonics,  we  obtain: 12 

1  00  r 

(p)  =  J2  j  ( lu> Wi  - (l + l)W>  tJ')dS  (61a> 

i  00  r 

=  ^^E(2/+  J  (WiUi)dS,  (61b) 

which  is  equivalent  to  the  formulae  (18)  in  Zschau  (1978)  and 
(5)  in  Platzman  (1984).  This,  however,  is  the  last  point  on  which 
we  are  still  in  agreement  with  our  predecessors. 

5.4.  Writing  the  Integral  as  a  Spectral  Sum 

Bringing  in  the  dynamical  Love  numbers  ki  and  the  phase 
lags  defined  in  Equation  (12),  one  can  express  the  products 
W/(t)  Ui(t)  via  the  spectral  components  of  the  disturbance 
W(r).13 

Although  the  formula 

OO 

£(2/ +  l)<Wi(0  *//(*)> 

1=2 

=  £(2/  +  l)f  W,\(0 )  km  Sin  e,(a>)  (62) 


where  =w  Wi(a>impq )  denotes  a  sum  of  all  terms  for  which 

co/mpq  takes  a  value  co. 

In  short:  first  sum  all  the  terms  corresponding  to  one  value 
of  co.  then  square  the  sum,  and  only  then  sum  over  all  the 
values  of  co. 

2.  Much  less  intuitive  is  the  fact  that  the  spectral  sum  will 
contain  extra  terms  that  are  missing  completely  in  the 
expression  (62).  As  we  shall  see  in  Appendix  F,  these 
terms  look  (up  to  some  caveat)  as  W/(co)  Wi(—oo).  They 
show  up  because  two  modes  of  opposite  values,  00  and  —co, 
correspond  to  the  same  physical  frequency  \a>\. 

For  the  time  being,  we  use  the  notation  ^  ' : 

OO 

£(2/+l  ){Wl{t)Ul{t)) 

1=2 

=  J^(21  +  !)  7  Wj2(<u)  ki(co)  sin  e,(&>),  (63) 

CO 

where  the  superscript  2  reminds  the  reader  that  the  spectral  sum 
needs  to  be  amended  down  the  road. 

Insertion  of  Equation  (63)  into  Equation  (61b)  results  in:15 


is  often  used  in  the  literature  (Zschau  1978;  Platzman  1984; 
Segatz  et  al.  1988), 14  accurate  examination  demonstrates  that 
it  is  incorrect.  To  appreciate  this,  one  simply  has  to  insert  the 
expansions  (2)  and  (12)  into  the  formula  (61b)  and  see  what 
happens. 

That  the  answer  differs  from  Equation  (62)  was  noticed  by 
Peale  &  Cassen  (1978).  However,  their  development  also  needs 
correction.  Later,  we  examine  this  matter  in  great  detail  and 
provide  a  full  inventory  of  the  terms  emerging  in  the  spectral 
expansion  for  damping  rate.  At  this  point,  we  only  mention  the 
two  key  circumstances: 

1.  The  conventional  expression  (62)  ignores  the  degeneracy 
of  modes,  that  is,  a  situation  where  several  modes  co/mpq 
with  different  sets  Impq  take  the  same  numerical  value  co. 
As  will  be  demonstrated  in  Section  6,  the  sum  over  modes 
co  in  Equation  (62)  should  be  substituted  with  a  sum  over 
distinct  values  of  the  modes: 


instead  of  Wfico)  in  Equation  (62)  use  this  : 

CO 


£ 


-,2 


^  ^  Wi((Oimpq ) 


U>lmpq=W 


12  On  the  boundary,  we  have:  Vl'(R)  =  EE2  WW  +  Ui(R)],  as  evident 
from  Equation  (59).  Together  with  Equation  (60),  this  expression  was  inserted 
in  Equation  (58).  By  doing  so,  we  omitted  the  diagonal  products  W\  Wi  and 
Ui  Ui  that  vanish  after  time  averaging.  (Indeed,  Wi  is  in  quadrature  with  W/, 
while  Ui  is  in  quadrature  with  f//.)  En  route  from  Equation  (61a)  to 
Equation  (61b),  we  took  into  account  that  the  time  averages  of  3 (t//  Wi)/dt 
also  vanish. 

13  In  this  subsection,  co  is  a  shortened  notation  for  the  mode  coimpq ,  not  the 
argument  of  the  pericenter. 

14  Our  expression  (62)  is  identical  to  the  upper  line  of  Equation  (10)  in 
Platzman  (1984).  Note  a  misprint  on  that  line  of  Platzman’s  equation:  a 
missing  factor  of  co.  Our  formula  (62),  when  truncated  to  /  =  2,  also  becomes 
equivalent  to  Equation  (22)  in  Zschau  (1978)  and  Equation  (12)  in  Segatz  et  al. 
(1988).  (Both  authors  kept  only  the  degree-2  terms.) 


(P)  =  ^^£(2/  +  l)®*/(©)sine/(fi>)  J  W,2(co)dS. 


(64) 


If  not  for  the  superscript  this  expression  would  coincide 
with  the  results  by  Zschau  (1978)  and  Platzman  (1984). 16  The 
superscript  reminds  us  of  the  important  caveat  in  the  evaluation 
of  the  sum:  the  factors  W2{co)  should  be  substituted  with  more 
complicated  expressions,  whereas  the  sum  should  be  carried  not 
over  all  modes  co  =  coimpq,  but  over  all  distinct  values  of  co,  see 
Section  6. 


6.  HEAT  PRODUCTION  OVER  TIDAL  MODES:  A 
SKETCHY  DERIVATION  OF  THE  FORMULA  (65),  IN 
NEGLECT  OF  THE  DEGENERACY 

We  must  insert  the  expansions  (2)  and  (12)  into  the  for¬ 
mula  (61b)  for  the  heating  rate,  in  order  to  obtain  a  comprehen¬ 
sive  version  of  the  somewhat  symbolic  sum  (Equation  (63))  and 
to  see  what  the  modified  sum  ^ 1  actually  means.  A  sketchy  ver¬ 
sion  of  this  calculation  (which  takes  into  account  that  the  modes 
may  have  either  sign,  but  neglects  the  degeneracy  of  modes)  is 
given  in  Appendix  F.  Extraordinarily  laborious,  the  full  calcula¬ 
tion  is  presented  in  Appendix  H  to  the  extended  version  of  this 
paper  (see  Efroimsky  &  Makarov  2014,  Appendix  H)  Here  we 
provide  the  final  results. 


15  Were  we  using  complex  potentials,  we  would  have  Wi  U*  instead  of  W;  (// 
in  Equation  (61b),  and  would  have  Wi  W*  instead  of  IV/  IV/  in  Equation  (63). 

16  Our  expression  (64)  should  be  compared  to  Equation  (22)  from  Zschau 

( 1978).  in  understanding  that  our  expression  furnishes  the  mean  damping  rate 
summed  over  the  entire  spectrum,  whereas  Zschau’s  formula  renders  the 
energy  loss  over  a  period  at  a  certain  frequency.  With  these  details  taken  into 
account,  the  formulae  are  equivalent.  They  are  also  equivalent  to  the  formulae 
(10)  and  (12)  in  Segatz  et  al.  (1988)  and  (10)  in  Platzman  (1984).  Note, 
however,  that  in  the  first  line  of  Platzman’s  formula  a  factor  of  co  is  missing. 
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In  the  case  of  a  uniformly  moving  pericenter,  the  average 
dissipation  rate  is: 


GM 


_„*2  ~  (R \2'+1  '  (/  -m)\  , 

w=— E(-)  <.(') 

1—2  v  7  m= 0  v  y  p= 0 


X  ^  Gfpqi.^  Xlmpq  kliXlmpq*)  sill(:i(Ximpq),  (65) 

q=—o o 


where  the  physical  frequencies  are  the  absolute  values  of  the 
Fourier  modes: 


Xltnpq  —  —  |(^  2p)  CO  +  (/  2 p  +  q)  j\cl 

+  m  (Q  —  0)|  ~  |(/  —  2 p  +  q)  n  —  md |,  (66) 

and  sine/(x/mp9)  is  what  they  often  call  1/2/  in  the  literature.17 

In  Efroimsky  &  Makarov  (2014,  Appendix  H),  we  have  also 
derived  a  formula  for  an  idle  pericenter,  but  the  applicability 
realm  of  that  formula  is  limited. 1 8 

Our  formula  (65)  differs  from  the  appropriate  expression  in 
Kaula  (1964,  Equation  (28))  that  contains  a  redundant  factor 
(l+*z)/2. 

In  the  special  situation  where 

1.  I  =  2, 

2.  the  body  is  incompressible,  so  k2  =  3  h 2/5, 19 

3.  the  spin  is  synchronous,  with  no  libration, 

the  expression  (65)  should  be  compared  to  the  formula  (31) 
from  Peale  &  Cassen  (1978).  This  comparison  is  carried  out 
in  Appendix  H.  In  our  expression,  all  terms  are  positive- 
definite,  because  the  factors  co2mpq  kii^impq)  sine2(®2mP9)  are 
even  functions  of  the  tidal  mode  co2mpq  ■  Peale  &  Cassen 
(1978),  however,  have  their  terms  proportional  to  the  products20 
(02 mpq  k2{(02mpq )  sin  |  e2(fi>2mpq )  |  which  are  negative  for  negative 
(Oimpq  ■  In  the  considered  setting,  the  largest  of  such  terms  were 
of  the  order  of  eA .  Such  inputs  lead  to  an  underestimation  of  the 
heat  production  rate  in  the  situations  where  the  eccentricity  is 
sufficiently  high  (like  in  the  case  of  the  Moon  whose  eccentricity 
might  attain  high  values  in  the  past  due  to  a  three-body  resonance 
with  the  Sun). 


7.  CONCLUSIONS 

We  have  derived  from  the  first  principles  a  formula  for  the 
tidal  dissipation  rate  in  a  homogeneous  spherical  body.  En  route 
to  that  formula,  we  compared  our  intermediate  results  with  those 
by  Zschau  (1978)  and  Platzman  (1984).  When  restricted  to  the 
special  case  of  an  incompressible  spherical  planet  spinning 
synchronously  without  libration,  our  final  formula  can  be 
compared  with  the  commonly  used  expression  from  Peale  & 
Cassen  (1978,  Equation  (31)).  The  two  turn  out  to  differ.  In 
our  expression,  the  contributions  from  all  Fourier  modes  are 
positive-definite,  which  is  not  the  case  of  the  formula  from 
Peale  &  Cassen.  As  a  result,  employment  of  the  formula  from 
Peale  &  Cassen  (1978)  may  cause  underestimation  of  tidal  heat 
production  in  some  situations.  For  example,  our  calculations 
for  the  rate  of  energy  dissipation  in  the  Moon  with  its  current 
parameters  yield  a  value  approximately  a  factor  of  two  greater 
than  the  result  derived  from  the  classic  equation,  with  the 
difference  coming  from  the  non-resonant  terms  having  the 
correct  sign. 

We  therefore  propose  to  use  our  Equation  (65)  for  rocky 
planets  and  moons,  instead  of  the  classic  formula  from  Peale 
&  Cassen  (1978,  Equation  (31)),  because  (1)  our  expression 
is  more  accurate  for  the  basic  case  of  objects  captured  into 
the  1:1  resonance,  and  (2)  it  correctly  captures  the  frequency 
dependence  of  tidal  dissipation  for  objects  outside  the  1:1 
resonance. 

Several  applications  are  provided  in  the  work  by  Paper  II, 
which  is  being  published  back  to  back  with  the  current  paper. 
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pointing  out  the  advantage  of  the  Lagrange  description,  and  to 
Gabriel  Tobie  for  a  very  useful  exchange  on  tidal  heating. 

The  authors  are  grateful  to  James  G.  Williams  for  a  meticulous 
reading  of  the  manuscript  and  very  judicious  comments  that 
were  of  great  help. 

The  authors’  special  thanks  go  to  the  referee,  Patrick  A. 
Taylor,  whose  thoughtful  and  comprehensive  report  enabled  the 
authors  to  improve  the  quality  of  the  paper  significantly. 
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17  It  would  not  hurt  to  reiterate  that  the  Fourier  modes  u>impq  can  be  of  either 
sign,  while  the  physical  forcing  frequencies  (Equation  (66))  are  positive 
definite.  Obviously,  Xlmpq  ki(ximpq)  sin  ei(ximpq)  =  a>lmpq  ki(a)impq) 

sin ei(a)[mpq),  because  the  dynamical  Love  numbers  are  even  functions, 
whereas  the  phase  lags  are  odd  and  of  the  same  sign  as  their  argument.  This  is 
why  the  tidal  quality  factors  may  be  expressed  as  1/(3/  =  sin  ^i(Xlmpq)  and 
also  as  1/(3/  =  |  sine i (coimpq)\,  with  the  absolute  value  symbols  being 
redundant  in  the  former  formula  and  needed  in  the  latter. 

18  For  an  idle  pericenter,  the  time-averaged,  tidal-heating  power  reads  as: 

{p)  =  agl  ^  ^  (2  _  Som )  E<p=o  F/mp(i) 

x  fZp'= 0  Pimp'd) 

ffq=— oo  P’lpq)e))P*lp'q')e)}q'=q—  2{p— p')  COs(2  (p  p)(P  o) 


APPENDIX  A 

THE  ASSOCIATED  LEGENDRE  FUNCTIONS  AND 
THEIR  NORMALIZATION 

The  Legendre  polynomials  are  usually  defined  by  the  Ro¬ 
driguez  formula: 

Id1,  , 

Pi{x)=-r-—l{x2-\)1.  (Al) 

2  l\  ax 

The  associated  Legendre  functions  Pim(x)  (termed  associated 
Legendre  polynomials  when  their  argument  is  sine  or  cosine  of 
some  angle)  were  introduced  by  Ferrers  (1877)  as21 


X  (Oimpq  k’li.Mlmpq')  sin  ^l((^lmpq)i 

co o  being  the  value  of  the  pericenter.  This  formula  is  of  a  limited  practical 
value,  since  coq  seldom  stays  idle.  For  example,  if  we  are  computing  tidal 
damping  in  a  planet  perturbed  by  the  star,  coq  of  the  star  as  seen  from  the  planet 
will  be  evolving  due  to  the  equinoctial  precession  of  the  planet  equator. 

19  Static  Love  numbers  of  an  incompressible  spherical  planet  satisfy  the 
relation  (21  +  1)  ki  =  3  hi.  As  explained  in  Appendix  G,  an  analogue  of  this 
equality  for  dynamical  Love  numbers  is  (21  +  1  )ki(coimpq)  =  3  hi(coimpq). 

20  In  the  notation  of  Peale  &  Cassen  (1978),  these  products  are  written  as 
(3/5)  h2  (2  -  2p  +  q  -  m)/Q2mPq- 


dm 

Plm(x)  =  (1  -  x2)m/2 - P,(x),  where  /  ^  m  f  0.  (A2) 

dxm 


21  Sometimes  in  the  literature  they  also  use  the  functions 
pf(x)  =  (-ir (i  -  x2r'2  p,(X)  =  (-i)m 

as  defined,  for  example,  in  Abramowitz  &  Stegun  (1972,  p.  332).  There  also 
exists  a  different  convention  wherein  P™(x)  lacks  the  (—  l)m  multiplier  and 
thus  coincides  with  Pim(x ),  as  in  Arfken  &  Weber  (1995,  p.  623). 
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The  so-defined  associated  Legendre  functions  are  sometimes 
called  unnormalized,  although  a  more  accurate  term  would  be 
in  Ferrers’  normalization.  This  normalization  reads  as: 

f1  2  0 l  +  m)\ 

/  Pim(x)Pim{x)dx  =  — - -  - - -Sw  (A3 a) 

y_i  21  +  1  (/  —  m)\ 

or,  equivalently: 


We  need  two  planetary  equations  in  the  Gauss  form  (Brouwer 
&  Clemence  1961,  p.  301,  Equation  (33)):22 


V 1  —  e2 


nae 


-R  cos  /  +  1  H - 

P 


sin(&>+  f)  coti  r 

— - — - W, 

nay/ 1  —  e2  a 


T  sin/ 


P,m  (sin  (p)Pr  (sin  /)  cos  (pdf 


2  (Z+m) ! 

21+  1  (/  —  m)\ 

(A3b) 


another  equivalent  form  being 


P/„,(C0S  (p)Pr,n( cos  <p)  sin  <pdcp 


2  (/  +  /»)! 

2Z+  1  (Z-m)! 

(A3c) 


The  associated  Legendre  functions  in  Ferrers’  normalization 
should  not  be  confused  with  the  associated  Legendre  functions 
P/m(x)  which  are  written  in  the  Schmidt  partial  normalization : 

J  Pim(x)Prm(x)dx  =  Yj-—^(2- 80m)8ir.  (A4) 

For  more  on  these  normalizations,  see  Winch  et  al.  (2005). 


APPENDIX  B 

KEEPING  cb  IMPLIES  EITHER  KEEPING  Mq 
OR  DEFINING  Ad  AS  dn/dt 

Under  disturbance,  the  mean  anomaly  is  written  as 

Ad  =  Mo  +  f  n(t)dt,  where  n(t )  =  yj G  (M  +  M*)a~2(t ) , 
Jt0 

(Bl) 


so  the  expression  (5)  for  the  Fourier  tidal  modes  acquires  the 
form  of 


tO/mpq  —  t'lllipq  ni  0  —  (Z  2 /?)  {(])  +  Ml q) 

+  q  Mq  +  (Z  —  2p  +  q)  n  +  m  (Q  —  9).  (B2) 

Kaula  (1964,  Equation  (40))  makes  an  oversight  by  accepting 
the  approximation 


Mo  =  -  ( cos  /  —  2  —  e  ^  A  —  (\+—  T  sin  /  , 

nae  LV  P  J  \  PJ 

where /is  the  true  anomaly,  p  =  a(l  —  e2)  is  the  semilatus 
rectum,  while  R,  T,  and  W  are  the  radial,  transversal,  and  normal- 
to-orbit  forces,  respectively.  In  a  situation  where  the  perturbation 
is  predominantly  transversal  and  the  terms  with  R  and  W  may 
be  neglected,  we  obtain: 

e  Vl  —  e2  (  r  \ 

co  +  Mo  » -  - -  (l+-  IT  sin/.  (B4) 

1  +  V 1  -  e2  na  \  p  ) 

A  low-inclined  moon  gets  predominantly  transversal  orbital 
disturbance  from  the  tides  it  creates  in  the  planet.  Inserting 
the  latter  expression  in  the  formula  (B2)  for  the  Fourier  modes, 
we  see  that  for  the  modes  with  a  zero  q  (like  the  semidiurnal 
tide  parameterized  with  Impq  =  2200)  the  input  from  the 
pericenter  rate  may  be  omitted  if  the  eccentricity  is  not  too 
large.  Indeed,  for  q  =  0,  the  term  qMo  vanishes,  while  the  term 
(Z  —  2p)(cb  +  Mo)  is  now  approximated  with  (Z  —  2p)  multiplied 
by  the  expression  (B4).  Although  co  and  Mo  can,  separately,  be 
substantial,  their  sum  (B4)  is  smaller  by  the  order  of  e.  Being 
(in  this  particular  case)  of  the  same  order  but  of  opposite  sign, 
co  and  Mo  largely  compensate  one  another.  Therefore,  if  we 
choose  to  drop  Mo,  we  should  also  drop  co.  In  this  special  case, 
dropping  both  will  be  legitimate. 

As  a  useful  aside,  we  would  remind  that  the  mean  longitude 
is  defined  through  L  =  M  +  co  +  £1,  its  rate  being  L  = 
y/G(M+  M*)a~i(t)  +  Mo  +  ch  +  Q.  As  we  have  just  seen, 
the  rates  Mo  and  co  largely  compensate  one  another  and  may 
both  be  neglected  in  the  considered  case.  If,  above  that,  the  rate 
of  the  node  happens  to  be  negligible,  then  the  mean  motion  from 
the  Kepler  law  will  be  close  to  the  mean  longitude  rate. 

B.2.  Example  2.  Orbital  Perturbation  Due  to  Oblateness 

The  situation  is  different  where  the  principal  perturbation  is 
due  to  the  oblateness  of  the  tidally  perturbed  primary.  The  mean 
rates  (Vallado  2007,  pp.  647-648) 


M  =  Ado  +n  n.  (B3) 

Indeed,  as  co  and  Mo  are  often  of  the  same  order,  it  is  incorrect 
to  keep  the  former  while  neglecting  the  latter.  We  present  two 
examples.  In  the  first,  co  and  Mo  are  of  the  same  order  but 
of  opposite  signs,  so  they  largely  compensate  one  another.  This 
suggests  a  simultaneous  neglect  of  both  terms.  In  the  second 
example,  co  and  Mo  turn  out  to  be  of  the  same  order  and  the 
same  sign,  so  keeping  one  of  these  terms  requires  keeping  the 
other. 

B.l.  Example  1.  Tidal  Perturbation  of  a  Low-inclination, 
Low-eccentricity  Orbit 

Consider  a  low-inclined  perturber.  From  the  tides  it  creates, 
the  perturber  gets  predominantly  transversal  orbital  disturbance. 


and 


3  R2  2  —  3  sin2  i 

Mo  =  —nJp  —  — - p— — 

4  a1  (1  —  e-yl2 


3  R2  4  —  5  sin2  i 
co  =  — nJ 2- 


a2  (1  —  e2) 


2\2 


are  of  the  same  order  and  sign.  Therefore,  when  keeping  co,  we 
must  also  include  Ado- 

The  easiest  way  to  get  rid  of  Ado  is  to  define  the  mean  motion 
as  h  =  Ad.  This,  the  so-called  anomalistic  mean  motion  will, 
however,  differ  from  yj G  ( M  +  M*)a~3(t) . 


22  The  system  (33)  in  Brouwer  &  Clemence  (1961,  p.  301)  contains  an 
equation  for  the  rate  de/dt,  where  (as  explained  on  the  preceding  page  in 
Brouwer  &  Clemence)  €  is  understood  as  € 1  =  A4o  +  co  =  A4q  +  co+  Q. 
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APPENDIX  C 

UNIVERSALITY  OF  THE  DARWIN-KAULA 
DESCRIPTION 

As  we  saw  earlier,  to  obtain  the  decomposition 
(12)  from  the  Fourier  series  (2),  each  term  of  the  latter  series 
must  be  endowed  with  a  mitigating  factor  k /  =  ki(a>[mpq)  of  its 
own  and,  likewise,  must  acquire  its  own  phase  lag  e;  =  €i(coimpq). 
In  the  literature,  some  authors  enquired  whether  this  mitigate- 
and-lag  method  is  general  enough  to  describe  tides.  The  answer 
to  this  question  is  affirmative,  insofar  as  the  tides  are  linear. 
Without  going  into  details  (to  be  found  in  Efroimsky  2012a, 
2012b),  we  would  mention  that  an  /-degree  part  of  the  operator 
(Equation  (11))  is  a  convolution  called  the  Love  operator. 

(R \  ^ 

-J  J  iq(t  r*,  t’)dt',  (Cl) 

Indeed,  linearity  of  tides  means  that,  at  each  time  t,  the  overall 
magnitude  of  reaction  depends  linearly  on  the  magnitudes  of 
the  disturbance  at  all  preceding  instants  of  time,  t'  ^  t.  The 
emergence  of  inputs  from  earlier  times  stems  from  the  inertia 
(“memory”)  of  the  material.  A  disturbance  that  took  place  at 
an  instant  t'  appears  in  the  integral  for  Ufr .  t)  with  a  weight 
left  —  t')  that  depends  on  the  elapsed  time.  Following  Churkin 
(1998),  who  gave  this  formalism  its  present  shape,  we  call  these 
weights  Love  functions. 

In  the  frequency  domain,  the  convolution  becomes: 

LR\l+l- 

U,(o>)  =  l-  J  k,{o>)  W,(co),  (C2) 

where  w  —  a>impq  is  the  tidal  mode  (not  the  periapse);  Ufa> )  and 
Wfco)  are  the  Fourier  images  or  Laplace  images  of  the  potentials 
Ufr ,  t)  and  W/(R,  r* ,  f);  while  the  complex  Love  numbers 

kfco)  =  \k,(co)\  e~ie‘<a,)  =  k,(co)  e~i(‘(w)  (C3) 

are  the  Fourier  or  Laplace  components  of  the  Love  functions 
left  —  t').  The  actual  dynamical  Love  numbers  are  the  real  parts 
of  the  complex  Love  numbers,  kfco)  =  |fc/(<y)|;  while  the  tidal 
lags  are  the  complex  Love  numbers’  negative  phases. 

The  frequency  dependencies  kfco)  and,  consequently,  kfco) 
and  6/  (cu)  can  be  derived  from  the  expression  for  the  complex 
compliance  J  (y)  or  the  complex  rigidity  fi  ( y )  =  I  /  J(y)  of  the 
mantle  (with  y  =  Ximpq  =  \(f>impq\  being  the  physical  forcing 
frequency).  The  dependency  J(x)  follows  from  the  rheological 
model. 

Evidently,  the  formula  (C2)  is  a  concise  version  of 
Equation  (12).  Thus  we  see  that  the  mitigate-and-lag  method 
ensues  directly  from  the  linearity  assumption. 

APPENDIX  D 

THE  EULERIAN  AND  LAGRANGIAN  DESCRIPTIONS. 
PERTURBATIVE  APPROACH  TO  A  PERIODICALLY 
DEFORMED  BODY 

D.l.  Perturbative  Treatment 

Under  perturbation,  two  changes  will  happen  in  a  point  r  at 
a  time  f. 

1 .  Physical  fields  will  now  acquire  different  values  in  this  point 
at  this  time. 


For  example,  a  moon  fixed  in  a  certain  position  relative 
to  the  planetary  surface  will  render  some  distribution  of 
its  potential  over  the  volume  of  the  host  planet.  The  same 
moon  fixed  in  a  different  position  will  generate  a  different 
distribution  of  its  potential.  This,  in  its  turn,  will  yield  a 
different  deformation  of  the  planet  and  therefore  a  different 
spatial  distribution  of  its  tidal-response  potential  and  of  all 
other  quantities. 

Thus,  instead  of  the  unperturbed  Eulerian  dependencies 
q°E(r ,  f),  we  now  have 

qE(r,t)  =  qE(r,  t)  +  q’E(r,t).  (Dl) 

2.  A  different  particle  will  now  arrive  in  the  point  r  at  the  time 
t.  It  will  not  be  the  same  particle  as  the  one  expected  there 
at  the  time  t  in  the  absence  of  perturbation. 

On  the  other  hand,  a  particle  that  starts  in  X  at  t  =  0,  will 
appear,  at  the  time  t,  not  in  the  point  x  —  f(X,  t)  but  in 
some  other  place 

r  =  x+  u  =  /(X,  t)  +  u(X,  t).  (D2) 

These  two  changes  will  influence  the  Lagrangian  dependencies 
on  the  initial  conditions.  The  dependency  of  each  field  will 
obtain  a  variation  q'L(X ,  t): 

qL(X,t)  =  q0L(X,t)  +  q'L(X,t).  (D3) 

In  the  absence  of  perturbation,  the  particle  X  was  destined  to  ar¬ 
rive  in  x,  wherefore  qE(X,  1)  was  defined  through  Equation  (17). 
Under  perturbation,  the  same  particle  X  is  expected  to  end  up 
in  r,  so  the  Lagrangian  dependency  becomes 

qL(X.  t )  =  qE(r,  t )  (D4a) 


=  qE(f(X ,  t)  +  u(X,  t ),  t)  (D4b) 


=  qE(f(X,  t),  t)  +  u(X,  t)Vx  qE  +  0( u2)  (D4c) 


=  qE(f(X,  t ),  t )  +  q'E(f(X,  t ),  t )  +  u(X,  f)V.v  qE  +  0( u2). 

(D4d) 

Insertion  of  Equation  (D3)  into  the  left-hand  side  of 
Equation  (D4)  will  give  us: 

q°L(X,  t)  +  q'L(X,  t )  =  qE(f(X.  t ),  t )  +  q'E(f(X,  t ),  t) 

+  u(X,  t)XxqE  +  <9(u2). 

Subtracting  Equation  (17b)  from  this  formula,  we  arrive  at 
a  relation  between  the  perturbations  of  the  Lagrangian  and 
Eulerian  quantities: 

q'L{X ,  t )  =  q'E(f{X,  t),  t)  +  u(X,  t)Vx  qE  +  0( u2),  (D5) 

where  the  first  term  on  the  right-hand  side,  qE(x,  t),  expresses 
the  change  in  the  final  spatial  distribution  of  the  field  q.  The  other 
two  terms  show  up  because  perturbation  changes  the  mapping 
from  X  to  the  current  location. 
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D.2.  An  Equivalent  Description 

A  slightly  different,  although  equally  valid  viewpoint  is 
possible.  In  a  reference  setting  at  time  t,  an  observer  located 
in  x  will  see  the  arrival  of  a  particle  that  started  from  X : 

q°L(X,  t )  =  q°E(x,  t).  (D6) 

In  a  perturbed  situation,  the  same  observer  in  jc  will  register,  at 
the  time  t ,  the  arrival  of  a  different  particle,  one  that  started  from 
X  -U: 

qL{X  -  U,  t)  =  qE{x,  t),  (D7a) 

which  is: 

qL(X ,  0  -  U VxqL  +  0(U2)  =  qE(x,  t).  (D7b) 

Subtraction  of  Equation  (D6)  from  Equation  (D7b)  gives  us 
the  variations: 

qL(X,  t )  -  q°L(X ,  t )  -  U VxqL  +  0(U2)  =  qE(x,  t )  -  qE{x,  t) 

(D8a) 

or,  simply, 

q'L(X,  t )  =  q'E(x,  t )  +  U VxqL  +  0(U2).  (D8b) 

Introducing  the  Jacobian  J  =  dV‘ /dV°  =  det(dxi/dXj), 
we  write: 

u  VxqL  =  U7  VA  qE  =  uV,  qE,  (D9) 

with  u  =  UJ.  Thus  Equation  (D8b)  and  Equation  (D5)  are 
equivalent  insofar  as  0( U2)  =  0(u2). 

While  the  language  of  Equation  (D5)  is  more  conventional 
than  that  of  Equation  (D8),  the  latter  description  is  easier  for 
physical  interpretation.  Suppose  we  are  observing  the  gradual 
cooling  of  a  flow.  In  an  unperturbed  setting,  a  particle  that 
started  in  X  at  the  time  t  —  0,  will  show  up  in  x  at  the  time  t. 
Accordingly,  a  measurement  of  the  temperature  in  x  at  the  time 
t  will  render,  in  the  absence  of  perturbation,  a  value  to  which 
the  particle  X  has  cooled  down  by  this  time — see  the  equality 
(Equation  (D6)). 

Under  perturbation,  the  rate  of  cooling  of  each  particle  will 
change.  In  addition,  owing  to  the  change  of  trajectories,  a 
different  particle  will  show  up  in  x  at  the  time  t.  Now  this 
will  be  a  particle  that  started  its  movement  at  t  —  0  from  some 
point  X  —  U.  So  a  measurement  of  the  temperature  in  the  point 
x  at  the  time  t  will  now  give  us  a  temperature  value  to  which  the 
particle  X  —  U  has  cooled  down — see  Equation  (D7a). 

The  difference  between  the  measurements  performed  in  x 
in  the  perturbed  and  unperturbed  cases  will,  according  to 
Equation  (D8b),  read  as  qrE(x,  t )  =  qE(X,  t)  —  \]VxqL  +  0(JJ2). 
The  first  term  on  the  right  renders  the  cooling  down  of  the 
particle  arriving  in  x  at  the  time  t,  while  the  second  and  third 
terms  reflect  the  fact  that,  under  disturbance,  we  register  a 
particle  arriving  from  a  point  displaced  by  U,  compared  to  the 
particle  that  would  be  brought  to  x  by  an  unperturbed  flow. 

With  aid  of  Equation  (D9),  the  expression  (D8b)  can  be 
equivalently  rewritten  as 

D,  =  d,+\VX,  (DIO) 

where  v  =  3u/3r,  while  D  =  d/dt  is  the  comoving  derivative. 
The  physical  interpretation  of  Equation  (DIO)  is  obvious:  the 


rate  of  cooling  of  a  moving  particle,  dq/dt,  can  be  measured  by 
a  quiescent  observer.  The  observer,  however,  must  amend  his 
result,  dq/dt,  with  a  correction  taking  into  account  the  fact  that, 
being  quiescent,  he  is  measuring  the  difference  between  the 
temperature  of  different  particles  passing  by,  not  of  the  same 
particle. 

D.3.  Periodically  Deformed  Solids.  Linearization 

Hereafter,  we  shall  restrict  our  consideration  to  the  case 
of  a  periodically  deformed  solid.  It  is  natural  to  associate 
the  reference  trajectory  x  =  f(X,  t)  with  the  equilibrium 
configuration.  In  this  configuration,  the  particles  stay  idle,  so 
x  coincides  with  the  initial  value  X : 

x  =  f(X,t),  where  f(X,  t)  =  X  for  all  t,  (Dll) 

while  all  properties  keep  in  time  their  fiducial  values: 

q0L(X,t)  =  q°E(x,t)  =  q0.  (D12) 

Be  mindful  that  the  superscript  0  did  not  originally  mark  a  value 
fixed  in  time,  but  a  trajectory  chosen  to  be  reference.  It  is  only 
now  that  the  role  of  a  reference  configuration  is  played  by  an 
equilibrium  body,  that  the  superscript  0  begins  to  denote  an 
unchanging  value. 

For  a  particle  originally  located  in  X,  its  perturbed  trajectory 
/•  differs  from  its  reference  trajectory  jc  by  some  u: 

r  =  x+u=  f(X,t)  +  u(X,t).  (D13) 

When  the  reference  trajectory  is  the  equilibrium,  insertion  of 
Equation  (Dll)  into  Equation  (D13)  results  in 

r  =  X  +  u(X,  t),  (D14a) 

which  can  also  be  written  as 

r  —  x  +  u(x,  t),  (D14b) 

because,  in  this  case,  the  unperturbed  trajectory  x  always 
coincides  with  the  initial  value  X. 

We  work  in  a  linearized  approximation,  neglecting  the  term 
<9(u2)  in  Equation  (D5)  and  writing  all  expansions  up  to  terms 
linear  in  the  displacement  u  or  velocity  v  =  du/dt. 

For  a  short  and  simple  explanation  of  the  linearized  La- 
grangian  and  Eulerian  descriptions  of  tides,  see  Wang  (1997).  A 
more  comprehensive  treatment  is  offered  in  the  book  by  Dahlen 
&  Tromp  (1998,  Section  3.1.1). 

D.4.  Conservation  of  Mass  in  the  Lagrangian 
and  Eulerian  Descriptions 

Denote  the  Eulerian  value  of  the  mass  density  with  pE(r,  t). 
As  mass  cannot  be  destroyed  or  created,  so  its  amount  in  a 
comoving  volume  V  of  a  flow  stays  constant: 

^  [  pE  dV‘  =  0.  (D15) 

dt  Jyt 

For  the  reference  history  x  =  f(X,  t),  this  would  imply: 

I  pE{x,t)d4x=[  pE(X,0)d4X.  (D16a) 

Jv;cf  Jv  o 
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For  a  perturbed  history  r  =  f(X,  t )  +  u(Y,  t),  we  have: 

f  pE(r,t)d4r  =  f  pE(X,0)d4X.  (D16b) 

Jv^  Jv° 

For  each  individual  particle,  its  perturbed  trajectory  r  stems 
from  the  same  initial  position  X  as  the  appropriate  reference 
trajectory  x,  so  the  initial  densities  are  the  same: 

Pe(X,  0)  =  pE(X,  0).  (D17) 

At  later  times,  however,  pE(r,t)  and  p^(x,t)  have  different 
functional  forms. 

D.4.1.  The  Continuity  Law  in  the  Eulerian  Description 

The  right-hand  sides  of  the  formulae  (D16a)  and  (D16b) 
coincide,  as  they  render  the  mass  of  the  same  initial  distribution 
pE(X,  0)  =  pE(X,  0).  Thus  the  left-hand  sides  of  the  two 
formulae  also  coincide: 

0=1  pE(r,t)d4r  —  j  pE(x,  t)d4x,  (D18) 

•'Vp'e*  JKt 

as  the  mass  stays  unchanged,  no  matter  whether  the  system 
follows  the  reference  history  or  a  perturbed  one.  Now  switch 
from  the  perturbed  coordinates,  r,  to  the  reference  ones,  jr: 

0  =  J  d4x[JpE{r ,  t )  -  p°(x,  0]  (D19a) 

=  J  t/4*[(l  +  Vt  •  u )pE(r,  t)  -  p°(x,  t)],  (D19b) 

where  the  Jacobian  is: 

dV‘  rt  dr,  9 

/  =  =  det—  =  1  +  Vr  •  u  +  0( u2) 

dKf  dxJ 

=  1  +  /V,.  -u+6»(u2).  (D20a) 

From  this,  we  see  that  the  Jacobian  can  also  be  written  as21 

1  +  0(u2)  , 

J  = - —  =  1  +  Vr  •  u  +  0(u2).  (D20b) 

1  —  V,  •  u 

From  Equation  (D19a),  we  obtain  the  exact  equality 

pE{f ,  t)J  =  pE(x,  t),  (D21) 

a  linearized  version  thereof  being 

0  =  pE(r ,  f)(l  +  V,.  •  u)  —  p£(r  —  u,  t)  +  0(u2)  (D22a) 

=  Pe(t,  t )  -  pg(r,  t)  +  pE{r,  t)Vr  ■  u  +  uVr  p%(r,  t)  +  0(u2). 

(D22b) 

For  a  small  t ,  the  deviation  u  between  the  two  trajectories  is 
linear  in  time,  and  so  is  the  difference  between  the  perturbed  and 
reference  density  functions.  Thus,  we  may  change  t)  Vr  u 

23  The  expression  (D20b)  indicates  that,  within  a  linear  approximation,  we  do 
not  need  to  distinguish  between  Vr  and  V*  =  /  Vr  =  Vr  +  (Vr  •  u)  Vr  + 
0(u2),  when  the  operators  are  applied  to  a  quantity  that,  by  itself,  is  of  the  first 
order  of  smallness. 
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to  p%(r,  t)  Vr  ■  u,  to  obtain  an  expression  correct  to  first  order 
in  u: 

P£  +  Vr- (p^u)  =0,  (D23) 

where  the  finite  variation  is 

Pe'(j,  t)  =  pE(r ,  t )  -  pE(r,  t).  (D24) 

We  would  reiterate  that  the  perturbative  approach  to  Eulerian 
quantities  implies  a  comparison  between  their  present  spatial 
distributions.  So  the  two  histories  are  compared  in  the  same 
point  r  and  at  the  same  time  t. 

In  Equation  (D22b),  we  could  also  have  changed  uVr  pE(r,  t ) 
to  uVr  pE(r ,  t ).  Then,  instead  of  Equation  (D23),  we  would  have 
obtained  pE  +  Vr  •  (pE  u)  =  0,  without  the  superscript  0  in  the 
second  term.  In  the  linear  approximation,  however,  this  would 
be  no  better  than  Equation  (D23).  Traditionally,  the  form  (D23) 
is  preferred  in  the  literature. 

However,  when  switching  to  a  differential  form  of  the 
conservation  law,  we  no  longer  need  to  keep  the  superscript  0, 
because  the  difference  between  pE(r ,  t)  and  pE(r,  t)  becomes 
infinitesimally  small.  So  the  differential  law  reads  as: 

+  V,-  ■  (pE  v)  =  0,  (D25) 

at 

where  v  =  3u/3 1,  the  partial  derivative  giving  the  rate  of  change 
with  coordinates  fixed. 

Employment  of  the  perturbative  formula  (D23)  near  a  de¬ 
formable  free  boundary  requires  some  care.  On  the  one  hand, 
the  reference  density  pE(r)  makes  an  abrupt  step  there.  On  the 
other  hand,  due  to  deformation  of  the  boundary,  we  may  get  a 
finite  present  density  in  a  point  where  the  reference  density  used 
to  be  zero,  and  vice  versa. 

D.4.2.  The  Continuity  Law  in  the  Lagrangian  Description 

The  Lagrangian  density  is  introduced  in  the  standard  way 
(Equation  (D4a)): 

pL{X,  t )  =  pE(r,  t),  (D26) 

so  the  formula  (D21)  becomes: 

pL{X,t)J  =  p“(x,0-  (D27a) 

For  a  periodically  deformed  solid,  the  reference  density  pE(x,  t) 
is  the  density  of  the  undeformed,  stable  configuration.  So 
Pg(x,t)  —  pE(X,  0)  =  p°(X)  is  time-independent,  and  the 
equality  (Equation  (D27a))  becomes  simply 

pL(X,t)J  =  p°E{X).  (D27b) 

In  accordance  with  the  general  formula  (D5),  we  interrelate 
the  density  variations  as 

p'L{X,  t)  —  p'E(x,  t )  +  u Vxp£(x,  t).  (D28a) 

In  the  considered  case  of  small  periodic  variations,  the  reference 
trajectory  is  simply  x  =  X  at  all  times;  so  on  the  right-hand 
side  of  the  above  formula  we  have  a  gradient  of  a  constant-in- 
time  stationary  distribution:  Vxp°( x,  t)  =  Vx p°(X).  Thence 
we  obtain: 

p'L(X,  t )  =  pE'{X ,  t )  +  uVx p°{X).  (D28b) 
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Combining  this  formula  with  Equation  (D23),  we  arrive  at24 


Pl  +  •  u  =  0, 

(D29) 

where  p'L  =  p'L(X,  t ),  while  pE  =  pE(X). 

D.5.  The  Poisson  Equation 

D.5.1.  In  the  Eulerian  Description 

Perturbed  or  not,  the  density  always  obeys  the  Poisson 
equation,  while  the  perturbing  potential  W  obeys  the  Laplace 
equation  outside  the  perturber: 

V,2  VE  =  -4  tzGpe, 

(D30a) 

V>r"  =  -47Gp“, 

(D30b) 

£ 

ii 

p 

(D30c) 

APPENDIX  E 

BOUNDARY  CONDITIONS 

The  boundary  condition  on  the  total  Eulerian  potential 
VEuier  is  trivial.  To  avoid  infinite  forces,  the  potential  must  be 
continuous: 

(exterior)  (interior) 

v£  =  v£  .  (Ei) 

The  boundary  condition  on  the  potential’s  gradient  emerges  as 
a  corollary  of  the  Gauss  theorem  and  therefore  mimics  a  similar 
condition  from  electrostatics.25  Let  a  small  area  s  =  s  n  of  the 
free  surface  be  sandwiched  between  the  top  and  bottom  of  a 
cylinder  of  an  infinitesimal  height  u  =  u  ■  s/s,  with  the  vector 
u  being  the  tidal  displacement.  The  top  and  bottom  should  each 
have  the  principal  curvature  radii  coinciding  with  those  of  the 
free  surface,  but  in  the  leading  order  this  can  be  ignored,  with 
the  enclosed  volume  thus  being  us  =  u  ■  s.  In  neglect  of  the 
contributions  from  the  infinitesimally  small  side  areas  of  the 
cylinder,  employment  of  the  Gauss  theorem  for  the  Eulerian 
potential  gives: 


Subtraction  of  Equation  (D30b)  from  Equation  (D30a)  results 
in  a  Poisson  equation  for  the  density  perturbation: 

V,.2  VE'  =  -4ttGpe'  (D31) 

or,  equivalently: 

V2  UE  =  -4jtGPe',  (D32) 

where  we  took  into  account  the  relations  (31)  and  (D30c). 


-4tt  Gp 


°“’=l 


=  /  Vx  VE  ■  ds 


=  vxy£ 


(exterior) 


■  s  -  vx  y£ 


(interior) 


•  s.  (E2) 

Over  a  surface  between  layers,  the  condition  will  read  as 


A  (exterior)  n  (interior) 

(4tt Gp°u  ■  s)  -(4itGp°u-s) 

(exterior)  (interior) 

=  Vx  VE  ■  S  -VXVE  ■  s  (E3) 

or,  equivalently, 


D.5.2.  In  the  Lagrangian  Description 

Insertion  of  the  formulae  (27)  and  (32)  into  the  Eulerian 
version  of  the  Poisson  equation,  Equation  (D31),  results  in  the 
Lagrangian  version  of  this  equation: 

Vr2  (y,/  -  u  •  Vx  V°)  =  -47 tG(p'l  -  u  •  Vx  p°).  (D33a) 

A  switch  to  differentiation  over  the  initial  position,  Vx,  would 
entail  corrections  of  the  order  of  0(u2).  In  neglect  of  those,  the 
equation  may  be  written  as 

V2(VL'  -  u  ■  Vx  V°)  =  —4tcG(p'l  -  u  •  Vx  p°).  (D33b) 


(exterior) 


-AitGp^u  h — -  VE 

3n 


-AitGp^u  h — -  VE 

3n 


(interior) 


(E4) 


In  application  to  tides,  it  can  be  interpreted  like  this:  the 
discontinuity  in  attraction  is  equal  to  the  attraction  of  the 
deformation  bulge  (Legros  et  al.  2006).  Since  yo,  W  and  their 
normal  gradients  are  continuous  on  the  boundary,  the  conditions 
on  U  and  V'  look  exactly  like  Equations  (E1)-(E4).  Specifically, 
in  Section  5.1  we  need  the  conditions  on  the  total  variation  V'\ 


For  an  initially  homogeneous  body,  the  above  formulae 
simplify  to: 

V,2(yL'  -  u  •  Vx  y  °)  =  —  4jt Gp'L  (D34a) 

and 

V2(  VL'  -  u  ■  Vx  y  °)  =  -4; r  Gp'h .  (D34b) 


(exterior)  (interior) 


(E5) 


(exterior) 


-47rGp°u  +  —  VE' 

3n 


(interior) 


-47rGp°u  +  —  VE' 

3n 


(E6) 


24  When  combining  Equation  (D28b)  with  Equation  (D23),  we  should  not  be 
confused  by  the  fact  that  in  Equation  (D28b)  all  quantities  are  functions  of  X, 
while  in  Equations  (D23)  and  (D24)  these  quantities  show  up  as  functions  of  r. 
Nor  should  we  be  confused  by  V  denoting  Vr  in  Equation  (D23)  and  Vx  in 
Equation  (D28b).  As  our  intention  is  simply  to  compare  the  functions,  we  can 
easily  change  the  notations  in  Equations  (D23)  and  (D24)  from  r  to  X, 
whereafter  Equation  (D29)  will  come  out  trivially. 


The  Eulerian  and  Lagrangian  potentials  are  interrelated  through 
VE'  =VL'-u-  Vx  y°.  (E7) 


25  Melchior  (1972)  attributes  the  derivation  of  the  boundary  condition  to 
Michel  Chasles. 
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Thence,  in  the  Lagrangian  description,  the  conditions  will 
acquire  the  form  of 


„  (exterior)  „ 

[VL'  -  u  ■  V,  V°]  =  [VL'  -  u  •  Vx  V°  ]  (E8) 

and 

(exterior) 


(interior) 


— 04' -u-  VJtV°)-4jrGp°u 

3n 


(interior) 


3n 


(Vl'-u- VXV°)  -  G  p°u 


(E9) 


When  the  boundary  is  welded  or  its  normal  is  parallel  to  VA  V  °, 
the  term  — u  ■  Vx  V  0  becomes  continuous  (Wang  1997).  It,  thus, 
can  be  removed  from  Equation  (E8),  rendering  the  incremental 
Lagrangian  potential  continuous.  This  term,  however,  cannot  be 
omitted  in  Equation  (E9). 


APPENDIX  F 

HEAT  PRODUCTION  AT  DIFFERENT  TIDAL  MODES. 

A  SKETCHY  DERIVATION  OF  THE  FORMULA  (65), 

IN  NEGLECT  OF  THE  “DEGENERACY” 

To  compute  the  dissipation  rate  at  separate  tidal  modes,  it 
is  necessary  to  insert  the  expansions  (2)  and  (12)  into  the  for¬ 
mula  (61b)  for  the  heating  rate.  This  will  render  a  comprehensive 
version  of  the  somewhat  symbolic  sum  (Equation  (63))  and  will 
enable  us  to  understand  what  the  modified  sum  :  actually 
means.  A  full  calculation  is  presented  in  Appendix  H  to  the 
extended  version  of  our  paper  (Efroimsky  &  Makarov  2014). 
Here  we  present  a  simplified  sketch  of  that  derivation. 

Recall  that  several  different  Fourier  modes  u>impq  can  share 
the  same  value  co.  Borrowing  a  term  from  quantum  mechanics, 
we  call  this  the  degeneracy  of  modes.  As  a  prelusory  exercise, 
we  calculate  dissipation  at  different  modes,  neglecting  the 
degeneracy.  In  other  words,  suppose  that  all  Fourier  modes  co  = 
oo\mpq  have  different  values.  Under  this  simplifying  assumption, 
the  expression  under  the  integral  in  Equation  (61b)  becomes: 

OO 

Y/(2l  +  D(W,(t)Ul(t)) 

1=2 

=  y  (2/  +  1  ){Wfco’)  cos [«' f  +  cpwfco’)]  (—ft>)  Ufa) ) 

(0,0)' 

x  sin[ft>?  +  (pufco)]) 

=  -  y<2/  +  1)|  Wfco')  t//(ft>)(sin[(ft)  -  co')t 

co,co' 

+  (Pufo)  —  (pwfa)')]  +  sin[(ft)  +  oo')t 
+  <Pu,(u)  +  <Pw,  (ft/)]},  (Fla) 

where  (•  ■  •)  denotes  time  averaging.  Of  the  two  sine  functions 
on  the  right-hand  side,  we  would  have  kept  only  the  first,  had  the 
Fourier  tidal  modes  been  positive-definite.  In  the  tidal  theory, 
however,  the  Fourier  modes  co  =  coimpq  can  assume  either  sign, 
so  both  sine  functions  must  be  taken  into  account: 

OO 

^(2/  +  l)(W,(0t>/(0> 

1=2 

ECO 

(2/  +  1)  -{  Wfco)  Ufo)  sin [(pufco)  -  (pwfco)] 

CO  ~ 

+  Wf—co)  Ufa))  sin [(pufco)  +  <pwf-co)]}  (Fib) 


=  (2/  +  1)  |  Wi2(°f  kfco)  sin  efa))  +  y  (21  +  1 )  | 

CO  CO 

x  Wfco)  Wf— co)  kfco)  sin  e'fco),  (Flc) 

where  we  recalled  that  the  dynamical  Love  number  is  an  even 
function  of  the  Fourier  mode. 

On  the  right-hand  side  of  Equation  (Flc),  the  first  sum  is 
an  expected  input  coinciding  with  the  expression  obtained  by 
other  authors — see,  for  example,  the  first  line  of  formula  (10)  in 
Platzman  (1984). 2(1  This  input  is  proportional  to  kfco)  sin  efco), 
where 

6/ (ft))  =  (pwfco)  -  (Pufco)  (F2) 

is  the  tidal  phase  lag  at  the  frequency  co  =  a)impq. 

The  second  sum  in  Equation  (Flc)  comes  into  being  due 
to  the  fact  that  the  Fourier  modes  are  not  positive-definite. 
This  input  contains  a  factor  of  kfco)  sin  e'fco),  where  the  angle 
e'fco)  is,  generally,  different  from  the  phase  lag  (Equation  (F2)) 
appropriate  to  the  mode  co  =  coimpq.  Indeed, 

effi))  =  —  [cpufco)  +  (pwf—o))]  =  (pwfco) 

—  (Puf(L>)  —  <Pwf«) )  —  I Pwf—co) 

=  efa))  -  [(Pwfw)  +  <pwf-u>)\.  (F3) 

At  first  glance,  this  result  is  most  unphysical.  Usually,  to  calcu¬ 
late  dissipation  rate,  we  have  to  sum,  over  physical  frequencies 
or  over  Fourier  modes,  terms  proportional  to  the  sines  of  phase 
lags  at  those  modes.  The  addition  of  a  finite  phase  to  those 
lags  looks  bizarre.  However,  an  accurate  calculation  carried  out 
in  Efroimsky  &  Makarov  (2014,  Appendix  H)  shows  that  the 
phase  consists  of  two  parts.  One  is  equal  to  [(— 1)/  —  1  |tt/2,  so 
its  presence  renders  an  overall  factor  of  (— 1);.  The  other  part  of 
the  phase  is  (in'  +  m)  X,  so  after  integration  over  the  surface,  it 
results  in  a  S(m ’  +  m)  factor,27  where  m  is  the  second  index  of 
coi,npq  =  co,  while  m'  is  the  second  index  of  coim'P’q <  =  —  co: 

i  00  r 

{p)  =  4 )Tgr^(2/+  1}  ] 

=  E<2/+  rif  /  dSWfco)[Wfco) 

CO 

+  (— 1  y<5(m'+  m)Wf— co)]kfco)  sine fco).  (F4) 

The  indices  m  and  in'  being  nonnegative  (see  Equation  (2)),  the 
emergence  of  S(m'  +  m)  indicates  that  the  summation  in  the 
second  part  must  be  reduced  to  m  =  in'  =  0: 

y,  —  (21  +  1  )k/  (co)  sin  e/  (co)  Wj2  (co) 

.CO=  COlmpq 

+  y  ~^(2l  +  l)kfa))  sine/(ft>)(— l/  Wfco)  Wf—co)  . 

CO=(Ol0pq 

(F5) 

26  The  second  line  in  Platzman’ s  formula  renders  oceanic  and  atmospheric 
inputs. 

27  The  finite  phase  assumes  the  value  of  [(—  1)*  —  1]  7r/2  +  {m'  +  m)k,  with 
the  integer  m'  being  the  order  of  a/  =  coim>p'q’,  and  m  being  that  of  <z>  =  coimpq. 
The  presence  of  [(—  1)*  —  1]  7r/2  in  the  phase  is  equivalent  to  multiplying  the 
sum  by  (—  1)*.  So,  for  even  /,  this  part  of  the  phase  can  be  ignored.  The 
presence  of  the  term  ( m '  +  m)  X  in  the  phase  tells  us  that,  after  integration  over 
the  surface,  only  the  terms  with  m'  =  m  =  0  stay.  Thus,  after  integration,  we 
are  effectively  left  with  €[(co)  =  (—  l/  €i(co)  8(m'  +  m). 


(P)  = 


i 


4  ttGR 


/ 


dS 
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We  then  see  what  the  superscript  ft  introduced  in  Equations  (63) 
and  (64)  actually  implies: 

O)  (J)=(J)lmpq 

+  J2  -(-1  )lW,(co)W,(-co),  (F6) 

(0=0)10  pq 


where  the  first  sum  on  the  right-hand  side  is  complete  (i.e.,  goes 
over  all  modes),  while  the  second  sum  is  only  over  the  modes 
with  a  vanishing  second  index. 

Now,  what  is  Wi(a>)7  Naively,  Wi(w)  =  Wi(a>impq)  should  be 
the  real  magnitude  of  the  term  Wimpq  of  Kaula’s  series  (2).  In 
reality,  we  have  degeneracy  of  modes,  so  in  each  of  the  two 
Fourier  series  (for  W  and  for  U)  we  first  must  group  together 
the  terms  corresponding  to  each  actual  value  of  mode,  and 
only  afterward  should  we  multiply  the  series  by  one  another 
and  perform  time  averaging.  This  calculation  is  presented  in 
Efroimsky  &  Makarov  (2014,  Appendix  H).  In  the  case  where 
the  apsidal  precession  is  uniform,  the  answer  is: 


(P) 


GM* 


T.(- 


1=2 


2/+1  / 


EE  O') 

m= 0 p= 0 


E2  \l  ill )  • 

G ,,  ,  |  (^  ^0 m)  tOimpqki(cOimpq) 

q=—oo  ' 

x  sin  €i(coimpq).  (F7) 


APPENDIX  G 

INTERRELATION  BETWEEN  DYNAMICAL  LOVE 
NUMBERS,  FOR  AN  INCOMPRESSIBLE 
HOMOGENEOUS  SPHERE 


For  an  incompressible  homogeneous  spherical  body,  the  static 
Love  numbers  read  as 


ki 


3  1 

2(1-  1)  1+  A, 


and  hi 


21+  1  1 
2(1-  1)  1  +  A,’ 


where 


A,  = 


(2 12+  4 /+  3  )p 

IgpR 


3  (2  Z2  +  41+  3) /r 
4/jrGp 2  R2 
3  (2 12  +  4 1  +  3) 
4/  7T  G  p2  R2  J 


(Gl) 


(G2) 


p  and  J  =  1/p  being  the  relaxed  rigidity  and  compliance,  and 
G  being  Newton's  gravity  constant.  The  formulae  (Gl)  yield  a 
well-known  relation  connecting  the  static  Love  numbers: 


(21+  l)ki  =  3hi.  (G3) 


Expressions  (Gl)  are  obtained  by  solving  a  system  comprising 
the  static  version  of  the  Second  Law  of  Newton  and  the  con¬ 
stitutive  equation  interconnecting  the  stress  and  strain  through 
the  rigidity  p.  A  wonderful  theorem,  called  the  correspondence 
principle  or  the  elastic-viscoelastic  analogy,  tells  us  that  in  many 
situations  the  dynamical  versions  of  the  Second  Law  of  Newton 
and  constitutive  equation,  when  written  in  the  frequency  domain 
as  algebraic  equations  for  operational  moduli,  mimic  the  static 
versions  of  these  equations.  In  order  for  this  correspondence  to 
take  place,  the  accelerations  and  inertial  forces  should  be  neg¬ 
ligibly  small  (see,  e.g.,  Appendix  B  to  Efroimsky  2012a).  In 


that  case,  the  complex  Love  numbers  ki(w)  and  hi(co)  will  be 
expressed  through  the  complex  operational  moduli  p  or  J  in  the 
same  algebraic  manner  as  the  static  k/  and  hi  are  expressed  via 
the  static  p  or  J .  Also  recall  that  the  static  expressions  (Gl)  were 
derived  under  an  extra  assumption  of  incompressibility.  If  this 
assumption  is  also  valid  in  the  dynamical  case,  then  the  complex 
ki(co)  and  hi(co)  are  expressed  through  the  complex  p  or  J  by 
formulae  mimicking  Equation  (Gl),  whence  an  expression  like 
Equation  (G3)  ensues  for  ki(w)  and  hi(co).  Its  imaginary  part 
will  read  as: 

(21+  l)ki(co)sine(a>)  —  3hi(co)sme(a>),  (G4) 

where  k,(m)  =  \ki(co)\,  hi(a>)  =  \h{(co)\  and  co  =  (ohnpq. 

To  draw  to  a  close,  we  would  emphasize  that  in  the  static 
expression  (G2)  the  letters  p  and  J  =  1/p  stand  for  the 
static  (relaxed)  values  of  the  rigidity  and  compliance.  In  a 
dynamical  analogue  of  this  expression,  the  same  letters  denote 
the  unrelaxed  values. 


APPENDIX  H 

COMPARISON  OF  OUR  RESULT  WITH  THAT 
OF  PEALE  &  CASSEN  (1978) 

It  would  be  instructive  to  compare  our  formula  (F7)  with  the 
classical  result  by  Peale  &  Cassen  (1978).  To  this  end,  three 
items  must  be  kept  in  mind. 

1.  Peale  &  Cassen  tacitly  assumed  that  averaging  should  be 
carried  out  not  only  over  the  tidal  period  but  also  over 
the  apsidal  period — this  can  be  understood  from  how  their 
formulae  (21)  transformed  into  Equation  (22).  This  is  why 
their  resulting  formula  (31)  is  appropriate  to  compare  with 
our  expression  (F7). 

2.  As  the  derivation  in  Peale  &  Cassen  was  intended  for  the 
incompressible  case  and  for  /  =  2  solely,  we  should  use, 
for  the  purpose  of  comparison,  the  equality  k2(a>2mPq )  = 
3h2(a>2mpq)/5  derived  in  Appendix  G. 

3.  In  Peale  &  Cassen,  only  the  case  of  synchronous  rotation 

was  addressed,  with  u>2mpq  —  (2  —  2p)co+(/2  —  2p  +  q)n  + 
m  (£2  —  0)  (2  —  2 p  +  q  —  m)  n,  where  n  is  the  apparent 

mean  motion  of  the  perturber.  Librations  were  ignored. 

Taking  all  this  into  account,  we  write,  for  the  purpose 
of  comparison,  an  appropriately  simplified  version  of  our 
expression  (F7): 


(synchronous)  /  p .  (incompress)  _ 

\rn= 2  — 


GM 


*2 


21+1  2  2 

t)  XX>L(i) 

m=0 p= 0 


X 


2  (2  -m)\  3 

/  ,  G2 pq(tt)  (  (2  )  0}2mpq~^lt2(t02mpq') 


q=— oo 

x  sin e2(co2mpq). 


(2  +  777)! 


(HI) 


As  the  time  lag  in  our  formula  (13a)  is  always  positive-definite, 
the  sign  of  the  phase  lag  €i(a>impq)  coincides  with  that  of  the 
tidal  mode  a>impq,  wherefore  the  product  a>impq  sinei(a)impq )  can 
always  be  written  down  as  a  product  of  absolute  values: 


C+jnptj  Sill  ( l{(V/fnpq)  —  Otlwpq  '  \  SU1  Cl(o.)fmpq)\  —  ,  (H2) 

(dlmpq 

where  1/Qimpq  =  |  sine;(®/mp9)|  is  the  inverse  quality  factor, 
and  ximpq  =  \u>impq\  is  the  positive-definite  physical  forcing 
frequency.  For  synchronous  spin  and  1  =  2,  the  forcing 


18 


The  Astrophysical  Journal,  795:6  (19pp),  2014  November  1 


Efroimsky  &  Makarov 


frequency  is  xi mpq  ^  1 2  —  2  p+ q  —  m  \  n ,  whence  the  quadrupole 
input  into  the  power  is: 


(synchronous)/  p\  (incompress)  _ 
\r>l= 2 


GM 


*2 


n  V  2/+1  2 

t)  £ 


n= 0 


,  vA  ,  (2  —  m)\ 

E^ipO')  E  (2  -  <W) 


p=0 


q=—oo 


(2  + hi)! 


3  |2  —  2/j  +  q  —  m\ 

x  —  «2 - H, 

5  Qlmpq 


(H3) 


with  an  absolute  value  in  the  numerator. 

Peale  &  Cassen  (1978)  had  in  their  formula  (31)  simply 
(2  —  2 p  +  q  —  m)  instead  of  |2  —  2  p  +  q  —  m\.  As  a  result, 
their  expression  for  dissipation  rate  contained  negative  inputs 
from  some  Fourier  modes  (i.e.,  for  some  sets  of  m  pq).  Being 
of  the  order  of  e4,  such  inputs  lead  to  an  underestimation  of 
the  heat  production  rate  in  situations  where  the  eccentricity  is 
high.  The  presence  of  such  inputs  in  formula  (31)  from  Peale  & 
Cassen  was  pointed  out  by  Makarov  (2013),  who  explored  tidal 
heat  production  in  the  Moon  in  the  cause  of  its  orbital  evolution. 
Presumably,  the  Moon  was  formed  much  closer  to  the  Earth  than 
it  is  today,  and  could  be  captured  into  a  three-body  resonance 
with  the  Sun,  driving  the  orbital  eccentricity  to  high  values  for 
a  limited  timespan  (Touma  &  Wisdom  1994). 
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